perm filename V2J.IN[1,DEK] blob sn#285513 filedate 1977-05-28 generic text, type C, neo UTF8
COMMENT ⊗   VALID 00017 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00003 00002	folio 415 galley 1
C00016 00003	folio 417 galley 2
C00031 00004	folio 419 galley 3
C00041 00005	folio 421 galley 4 WARNING: Much of this tape was unreadable!
C00049 00006	folio 423 galley 5
C00061 00007	folio 426 galley 6
C00076 00008	folio 429 galley 7
C00092 00009	folio 432 galley 8
C00106 00010	folio 434 galley 9
C00121 00011	folio 439 galley 10 WARNING: Some bad spots were on this tape.
C00133 00012	folio 441 galley 11
C00149 00013	folio 444 galley 12
C00164 00014	folio 448 galley 13
C00176 00015	folio 451 galley 14
C00188 00016	folio 460 galley 15
C00198 00017	folio 461 galley 16
C00222 ENDMK
C⊗;
folio 415 galley 1
    0  {U0}{H10L12M29}58320#Computer Programming!(A.-W./Knuth)!f. 
    2  415!ch. 4!g. 1|'{A12}|π|∨4|∨.|∨5|∨.|∨2|9|∨T|∨h|∨e|9|∨G|∨r|∨e
    5  |∨a|∨t|∨e|∨s|∨t|9|∨C|∨o|∨m|∨m|∨o|∨n|9|∨D|∨i|∨v|∨i|∨s|∨o|∨r|'
    6  {A6}If |εu |πand |εv |πare integers, not both 
   14  zero, we say that their |εgreatest common divisor, 
   22  |πgcd(|εu,|4v), |πis the largest integer which 
   28  evenly divides both |εu |πand |εv. |πThis de_nition 
   36  makes sense, because if |εu|4|=|↔6α=↓|40 |πthen 
   42  no integer greater than |¬G|εu|¬G |πcan evenly 
   49  divide |εu, |πbut the integer 1 does divide both 
   58  |εu |πand |εv; |πhence there must be a largest 
   67  integer that divides them both. When |εu |πand 
   75  |εv |πare both zero, every integer evenly divides 
   83  zero, so the above de_nition does not apply; 
   91  it is convenient to set|'|h|πgcd(|εu,|4v)|4|∂α=↓|4|πgcd(|→α_
   96  ↓|εu,|4v),|E|n|;{A6}| |πgcd(0,|40)|4|Lα=↓|40.|J!(1)>
   98  {A6}!|4|4|4|4The de_nitions just given obviously 
  103  imply that|'{A6}| gcd(|εu,|4v)|4|Lα=↓|4|πgcd(|εv,|4u),|J!(2)
  105  >{A6}| |πgcd(|εu,|4v)|4|Lα=↓|4|πgcd(|→α_↓|εu,|4v),|J!(3)>
  107  {A6}| |πgcd(|εu,|40)|4|Lα=↓|¬Gu|¬G.|J!(4)>{A6}|π!|4|4|4|4In 
  109  the previous section, we reduced the problem 
  116  of expressing a rational number in ``lowest terms'' 
  124  to the problem of _nding the greatest common 
  132  divisor of its numerator and denominator. Other 
  139  applications of the greatest common divisor have 
  146  been mentioned for example in Sections 3.2.1.2, 
  153  3.3.3, 4.3.2, 4.3.3. So the concept of the greatest 
  162  common divisor is important and worthy of serious 
  170  study.|'*?*?!|4|4|4|4The |εleast common multiple 
  175  |πof two integers |εu |πand |εv, |πwritten lcm(|εu,|4v), 
  183  |πis a related idea which is also important. 
  191  It is de_ned to be the smallest positive integer 
  200  which is a multiple of (i.e., evenly divisible 
  208  by) both |εu |πand |εv; |πand lem(0,|40)|4α=↓|40. 
  215  The classical method for teaching children how 
  222  to add fractions |εu/u|¬S|4α+↓|4v/v|¬S |πis to 
  228  train them to _nd the ``least common denominator,'' 
  236  which is lcm(|εu|¬S,|4v|¬S).|'|π!|4|4|4|4According 
  240  to the ``fundamental theorem of arithmetic'' 
  246  (proved in exercise 1.2.4<21), each positive 
  252  integer |εu |πcan be expressed in the form|'{A6}|εu|4α=↓|42|
  260  gu|r23|gu|r35|gu|r57|gu|r711|gu|r1|r1|4|¬O|4|¬O|4|¬O|4α=↓|4|
  260  ≥u|uc|)p|4|πprime|)|4|εp|gu|rp,|J!(5)|;{A6}|πwhere 
  262  the exponents |εu|β2, u|β3,|4.|4.|4. |πare uniquely 
  268  determined nonnegative integers, and where all 
  274  but a _nite number of the exponents are zero. 
  283  From this canonical factorization of a positive 
  290  integer, it is easy to discover one way to compute 
  300  the greatest common divisor of |εu |πand |εv: 
  308  |πBy (2), (3), and (4) we may assume that |εu 
  318  |πand |εv |πare positive integers, and if both 
  326  of them have been canonically factored into primes, 
  334  we have|'{A6}gcd(|εu,|4v)|4|∂α=↓|4|≥u|uc|)p|4|πprime|)|4|εp|
  336  ur|πmin(u|βp,v|βp)|)|),|J!(6)|;{A6}| lcm(|εu,|4v)|4|Lα=↓|4|≥
  337  u|uc|)p|4|πprime|)|4|εp|ur|πmax(|εu|βp,v|βp)|)|).|J!(7)>
  338  |πThus, for example, the greatest common divisor 
  345  of |εu|4α=↓|47000|4α=↓|42|g3|4|¬O|45|g3|4|¬O|47 
  347  |πand |εv|4α=↓|44400|4α=↓|42|g4|4|¬O|45|g2|4|¬O|411 
  349  |πis 2|urmin(3,4)|)|)5|urmin(3,2)|)|)7|urmin(1,0)|)|)11|urmi
  350  n(0,1)|)|)|4α=↓|42|g3|4|¬O|45|g2|4α=↓|4200. The 
  352  least common multiple of the same two numbers 
  360  is 2|g4|4|¬O|45|g3|4|¬O|47|4|¬O|411|4α=↓|4154000.|'
  362  !|4|4|4|4From formulas (6) and (7) we can easily 
  370  prove a number of basic identities concerning 
  377  the gcd and the lcm:|'{A6}|hlcm(gcd(|εu,|4v),|4|πgcd(|εu,|4w
  382  ))|4|∂α=↓|4|πgcd(|εu,|4v)|4.|4|πlcm(|εu,|4v),!!|∂|πif!!|∂|εu
  382  ,|4v|4|¬R|40;|J!(10)|E|n|'| |πgcd(|εu,|4v)w|4|Lα=↓|πgcd(|εuw
  383  ,|4vw),|L|πif|L|εw|4|¬R|40;|J!(8)>{A6}| |πlcm(|εu,|4v)w|4|Lα
  384  =↓|4|πlcm(|εuw,|4vw),|L|πif|L|εw|4|¬R|40;|J!(9)>
  385  {A6}| u|4|¬O|4v|4|Lα=↓|4|πgcd(|εu,|4v)|4|¬O|4|πlcm(|εu,|4v),
  385  |L|πif|L|εu,|4v|4|¬R|40;|J!(10)>{A6}| |πgcd{H12}({H10}lcm(|ε
  386  u,|4v),|4|πlcm(|εu,|4w){H12}){H10}|4|Lα=↓|4|πlcm{H12}({H10}u
  386  ,|4gcd(|εv,|4w){H12}){H10};|J!(11)>{A6}{A6}|πThe 
  388  latter two formulas are ``distinctive laws'' 
  394  analogous to the familiar identity |εuv|4α+↓|4uw|4α=↓|4u(v|4
  399  α+↓|4w). |πEquation (10) reduces the calculation 
  405  of gcd(|εu,|4v) |πto the calculation of lcm(|εu,|4v), 
  412  |πand conversely.|'{A12}|≡E|≡u|≡c|≡l|≡i|≡d|≡'|≡s 
  415  |≡a|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m|≡.|9|4Although Eq. 
  417  (6) is useful for theoretical purposes, it is 
  425  generally no help for calculating a greatest 
  432  common divisor in practice, because it requires 
  439  that we _rst determine the factorization of |εu 
  447  |πand |εv. |πThere is no known method for _nding 
  456  the prime factors of an integer very rapidly 
  464  (see Section 4.5.4). But fortunately there is 
  471  an e∃cient way to calculate the greates common 
  479  divisor of two integers without factoring them, 
  486  and, in fact, such a method was discovered over 
  495  2250 years ago; this is ``Euclid's algorithm,'' 
  502  which we have already examined in Sections 1.1 
  510  and 1.2.1.|'!|4|4|4|4Euclid's algorithm is found 
  516  in Book 7, Propositions 1 and 2 of his |εElements 
  526  (c. 300|4{H7}|πB{H10}.{H7}C{H10}.), but it probably 
  531  wasn't his own invention. Scholars believe that 
  538  the method was known up to 200 years earlier, 
  547  at least in its subtractive form, and it was 
  556  almost certainly known to Eudoxus  |εc. 375 {H7}|πB{H10}.{H7
  564  }C{H10}. [Cf. K. von Fritz, |εAnn. Math. (2) 
  572  |≡4|≡6 (1945), 242<264.] |πWe might call it the 
  580  granddaddy of all algorithms, because it is the 
  588  oldest nontrivial algorithm which has survived 
  594  to the present day. (The chief rival for this 
  603  honor is perhaps the ancient Egyptian method 
  610  for multiplication, which was based on doubling 
  617  and adding, and which forms the basis for e∃cient 
  626  calculation of |εn|πth powers as explained in 
  633  Section 4.6.3. But the Egyptian manuscripts merely 
  640  give examples which are not completely systematic, 
  647  and which were certainly not stated systematically; 
  654  the Egyptian method is therefore not quite deserving 
  662  of the name ``algorithm.'' Several ancient Babylonian 
  669  methods, for doing such things as solving certain 
  677  sets of quadratic equations in two variables, 
  684  are also known. Genuine algorithms are involved 
  691  in this case, not just special solutions to the 
  700  equations for certain input parameters; even 
  706  though the Babylonians invariably presented each 
  712  method in conjunction with an example worked 
  719  with special parameters, they regularly explained 
  725  the general procedure in the accompanying text. 
  732  [See D. E. Knuth, CACM |≡1|≡5 (1972), 671677.] 
  740  Many of these Babylonian algorithms predate Euclid 
  747  by 1500 years, and they are the earliest known 
  756  instances of written procedures for mathematics. 
  762  They do not have the stature of Euclid's algorithm, 
  771  however, since they do not involve iteration 
  778  and since they have been superseded by modern 
  786  algebraic methods.)|'!|4|4|4|4Since Euclid's 
  790  algorithm is there fore important for historical 
  797  as well as practical reasons, let us {U0}{H10L12M29}58320#Co
folio 417 galley 2
  804  mputer Programming!(A.-W./Knuth)!f. 417!ch. 4!g. 
  808  2|'{A12}{H9L11}{I1.7L}|≡P|≡r|≡o|≡p|≡o|≡s|≡i|≡t|≡i|≡o|≡n|≡.|9
  809  |4|εGiven two positive integers, ⊂nd their greatest 
  816  common divisor.|'{A12}|πLet |εA, C |πbe the two 
  824  given positive integers; it is required to _nd 
  832  their greatest common divisor. If |εC |πdivides 
  839  |εA, |πthen |εC |πis a common divisor of |εC 
  848  |πand |εA, |πsince it also divides itself. And 
  856  it clearly is in fact the greatest, since no 
  865  greater number than |εC |πwill divide |εC.|'{A12}|πBut 
  873  if |εC |πdoes not divide |εA, |πthen continually 
  881  subtract the lesser of the numbers |εA, C |πfrom 
  890  the greater, until some number is left which 
  898  divides the previous one. This will eventually 
  905  happen, for if unity is left, it will divide 
  914  the previous number.|'{A12}Now let |εE |πbe the 
  922  positive remainder of |εA |πdivided by |εC; |πlet 
  930  |εF |πbe the positive remainder of |εC |πdivided 
  938  by |εE; |πand let |εF |πbe a divisor of |εE. 
  948  |πSince |εF |πdivides |εE |πand |εE |πdivides 
  955  |εC|4α_↓|4F, F |πalso divides |εC|4α_↓|4F; |πbut 
  961  it also divides itself, so it divides |εC. |πAnd 
  970  |εC |πdivides |εA|4α_↓|4E; |πtherefore |εF |πalso 
  976  divides |εA|4α_↓|4E. |πBut it also divides |εE; 
  983  |πtherefore it divides |εA. |πHence it is a common 
  992  divisor of |εA |πand |εC.|'{A12}|πI now claim 
 1000  that it is also the greatest. For if |εF |πis 
 1010  not the greatest common divisor of |εA |πand 
 1018  |εC, |πsome larger number will divide them both. 
 1026  Let such a number be |εG.|'{A12}|πNow since |εG 
 1035  |πdivides |εC |πwhile |εC |πdivides |εA|4α_↓|4E, 
 1041  G |πdivides |εA|4α_↓|4E. G |πalso divides the 
 1048  whole of |εA, |πso it divides the remainder |εE. 
 1057  |πBut |εE |πdivides |εC|4α_↓|4F; |πtherefore 
 1062  |εG |πalso divides |εC|4α_↓|4F. G |πalso divides 
 1069  the whole of |εC, |πso it divides the remainder 
 1078  |εF; |πthat is, a greater number divides a smaller 
 1087  one. This is impossible.|'{A12}Therefore no number 
 1094  greater than |εF |πwill divide |εA |πand |εC, 
 1102  |πand so |εF |πis their greatest common divisor.|'
 1110  {A12}|≡C|≡o|≡r|≡o|≡l|≡l|≡a|≡r|≡y|≡.|9|4This argument 
 1112  makes it evident that any number dividing two 
 1120  numbers divides their greatest common divisor. 
 1126  |εQ.E.D.|'{A12}{IC}{H10L12}Note.|9|4|πEuclid's 
 1128  statements have been simpli_ed here in one nontrivial 
 1136  respect: Greek mathematicians did not regard 
 1142  unity as a ``divisor'' of another positive integer. 
 1150  Two positive integers were either both equal 
 1157  to uhity, or they were relatively prime, or they 
 1166  had a greawt*?*?*?hey were relatively prime, or 
 1173  they had a greatest common divisor. In fact, 
 1181  unity was not even considered to be a ``number,'' 
 1190  and zero was of course nonexistent. These rather 
 1198  awkward conventions made it necessary for Euclid 
 1205  to duplicate much of his discussion, and he gave 
 1214  two separate propositions which each are essentially 
 1221  like the one appearing here.|'!|4|4|4|4In his 
 1228  discussion, Euclid _rst suggests subtracting 
 1233  the smaller of the two current numbers from the 
 1242  larger repeatedly, until we get two numbers in 
 1250  which one is a multiple of another; but in the 
 1260  proof he really relies on taking the remainder 
 1268  of one number divided by another (and since he 
 1277  has no concept of zero, he cannot speak of the 
 1287  remainder when one number divides the other). 
 1294  So it is reasonable to say that he imagines each 
 1304  |εdivision |π(not the individual subtractions) 
 1309  as a single step of the algorithm, and hence 
 1318  an ``authentic'' rendition of his algorithm can 
 1325  be phrased as follows:|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m 
 1330  |≡E (|εOriginal Euclidean algorithm).!|πGiven 
 1334  two integers |εA and |εC |πgreater than unity, 
 1342  this algorithm _nds their greatest common divisor.|'
 1349  {I1.8H}|≡E|≡1|≡.|9[|εA |πdivisible by |εC?] |πIf 
 1354  |εC |πdivides |εA, |πthe algorithm terminates 
 1360  with |εC |πas the answer.|'{A3}|≡E|≡2|≡.|9[Replace 
 1366  |εA |πby remainder.] If |εA |πmod |εC |πis equal 
 1375  to unit, the given numbers were relatively prime, 
 1383  so the algorithm terminates. Otherwise replace 
 1389  the pair of values |ε(A,|4C) |πby |ε(C,|4A|4|πmod|4|εC) 
 1396  |πand return to step E1.|'{IC}{A12}!|4|4|4|4The 
 1402  ``proof'' Euclid gave, which is quoted above, 
 1409  is especially interesting because it is, of course, 
 1417  not really a proof at all*3 He veri_es the result 
 1427  of the algorithm only if step E1 is performed 
 1436  once or thrice. Surely he must have realized 
 1444  that step E1 could take place more than three 
 1453  times, although he made no mention of such a 
 1462  possibility. Not having the notion of a proof 
 1470  by mathematical induction, he could only give 
 1477  a proof for a _nite number of cases. (In fact, 
 1487  he often proved only the case |εn|4α=↓|43 |πof 
 1495  a theorem which he wanted to establish for general 
 1504  |εn.) |πAlthough Euclid is justly famous for 
 1511  the great advances he made in the art of logical 
 1521  deduction, techniques for giving valid proofs 
 1527  by indukc*?*?*?e hot discovered until many centuries 
 1534  later, and the crucial ideas for proving the 
 1542  validity of |εalgorithms |πare only now becoming 
 1549  really clear. (See Section 1.2.1 for a complete 
 1557  proof of Euclid's algorithm, together with a 
 1564  short discussion of general proof procedures 
 1570  for algorithms.)|'!|4|4|4|4It is worth noting 
 1576  that this algorithm for _nding the greatest common 
 1584  divisor was chosen by Euclid to be the very _rst 
 1594  step in his development of the theory of numbers. 
 1603  The same order of presentation is still in use 
 1612  today in modern textbooks. Euclid also gave (Proposition 
 1620  34) a method to _nd the least common multiple 
 1629  of two integers |εu |πand |εv, |πnamely to divide 
 1638  |εu |πby gcd(|εu,|4v) |πand to multiply by |εn; 
 1646  |πthis is equivalent to Eq. (10).|'!|4|4|4|4If 
 1653  we avoid Euclid's bias against the numbers 0 
 1661  and 1, we can reformulate Algorithm E in the 
 1670  following way:|'{A12}|≡A|≡l|≡g|≡o|≡r|≡o|≡t|≡h|≡m 
 1673  |≡A (|εModern Euclidean algorithm).!|πGiven nonnegative 
 1678  integers |εu |πand |εv, |πthis algorithm _nds 
 1685  their greatest common divisor. (|εNote|*/: |\|πThe 
 1691  greatest common divisor of |εarbitrary |πintegers 
 1697  |εu |πand |εv |πmay be obtained by applying this 
 1706  algorithm to |ε|¬Gu|¬G |πand |ε|¬Gv|¬G, |πbecause 
 1712  of Eqs. (2) and (3){H12})|'{A6}{H10}|≡A|≡1|≡.|9|ε[v|4α=↓|40?
 1717  ] |πIf |εv|4α=↓|40, |πthe algorithm terminates 
 1723  with |εu |πas the answer.|'{A3}{I1.8H}|≡A|≡2|≡.|9[Take 
 1729  |εu |πmod |εv.] |πSet |εr|4|¬L|4u |πmod |εv, 
 1736  u|4|¬L|4v, v|4|¬L|4r, |πand return to A1. The 
 1743  operations of this step decrease the value of 
 1751  |εv, |πbut leave gcd(|εu,|4v) |πunchanged.)|'
 1756  {A12}!|4|4|4|4For example, we may calculate gcd(40902, 
 1762  24140) as follows:|'{IC}{A12}|h|πgcd(88888,|488888)|4|∂α=↓|4
 1765  gcd(0000,|40000)|4α=↓|4gcd(9999,|40000)|4α=↓|4gcd(1111,|4111
 1765  )|E|n|;| gcd(40902,|424140)|4|Lα=↓|4gcd(24140,|416762)|4α=↓|
 1766  4gcd(16762,|47378)>{A4}|Lα=↓|4gcd(7378,|42006)|4α=↓|4gcd(200
 1767  6,|41360)|4α=↓|4gcd(1360,|4646)>{A4}|Lα=↓|4gcd(646,|468)|4α=
 1768  ↓|4gcd(68,|434)|4α=↓|4gcd(34,|40)|4α=↓|434.>{A6}!|4|4|4|4A 
 1770  proof that Algorithm A is valid follows readily 
 1778  from Eq. (4) and the fact that|'{A6}gcd(|εu,|4v)|4α=↓|4|πgcd
 1785  (|εv,|4u|4α_↓|4qv),|J!(13)|;{A6}|πif |εq |πis 
 1789  any integer. Equation (13) holds because any 
 1796  common divisor of |εu |πand |εv |πis a divisor 
 1805  of both |εv |πand |εu|4α_↓|4qv, |πand, conversely, 
 1812  any common divisor of |εv |πand |εu|4α_↓|4qv 
 1819  |πmust divide both |εu |πand |εv.|'!|4|4|4|4|πThe 
 1826  following |¬m|¬i|¬x program illustrates the fact 
 1832  that Algorithm A can easily be implemented on 
 1840  a computer:|'{A12}|≡P|≡r|≡o|≡g|≡r|≡a|≡m |≡A (|εEuclid's 
 1845  algorithm).!|πAssume that |εu |πand |εv |πare 
 1851  single-precision, nonnegative integers, stored 
 1855  respectively in locations |¬u and |¬v; this program 
 1863  puts gcd(|εu,|4v) |πinto rA.|'{A6}{H9L11M29}|∂!!|4|∂333|4|4|
 1867  ∂!!!|4|∂!!!!!|∂!!!!!!!!!|∂|E|;|π|>|;|¬l|¬d|¬x|'
 1871  |¬u|'1|;rX|4|¬L|4|εu.|'>|>|π|;|¬j|¬m|¬p|'|¬2|¬f|'
 1879  1|;|;>|>|¬1|¬h|'|¬s|¬t|¬x|'|¬v|'|εT|;v|4|¬L|4|πrX.|'
 1888  >|>|;|¬s|¬r|¬a|¬x|'|¬5|'|{U0}{H10L12M29}58320#Computer 
folio 419 galley 3
 1894  Programming!(A.-W./Knuth)!f. 419!ch. 4!g. 3|'
 1898  {A12}The running time for this program is 19|εT|4α+↓|46 
 1906  |πcycles, where |εT |πis the number of divisions 
 1914  performed. The discussion in Section 4.5.3 shows 
 1921  that we may take |εT|4α=↓|40.842766|4|πln |εN|4α+↓|4|π0.06 
 1927  |πas an approximate average value when |εu |πand 
 1935  |εv |πare independently and uniformly distributed 
 1941  in the range |ε1|4|¬E|4u,|4v|¬E|4N.|'{A12}|π|≡A 
 1946  |≡b|≡i|≡n|≡a|≡r|≡y |≡m|≡e|≡t|≡h|≡o|≡d|≡.|9|4Since 
 1948  Euclid's patriarchal algorithm has been used 
 1954  for so many centuries, it is a rather surprising 
 1963  fact that it may not always be the best method 
 1973  for _nding the greatest common divisor after 
 1980  all*3 A quite di=erent gcd algorithm, which is 
 1988  primarily suited to binary arithmetic, was discovered 
 1995  by J. Stein in 1961 [see |εJ. Comp. Phys. |≡1 
 2005  (1967), 397<405.]|'{A6}|π|4|1|1|1a)|9If |εu |πand 
 2010  |εv |πare both even, then gcd(|εu,|4v)|4α=↓|42|4|πgcd(|εu/2,
 2015  |4v/2). |π[See Eq. (8).]|'|4|1|1b)|9If |εu |πis 
 2022  even and |εv |πis odd, then gcd(|εu,|4v)|4α=↓|4|πgcd(|εu/2,|
 2028  4v). |π[See Eq. (6).]|'|4|4c)|9As in Euclid's 
 2035  algorithm, gcd(|εu,|4v)|4α=↓|4|πgcd(|εu|4α_↓|4v,|4v). 
 2037  |π[See Eqs. (13), (2).]|'|4|1|1d)|9If |εu |πand 
 2044  |εv |πare both odd, then |εu|4α_↓|4v |πis even, 
 2052  and |ε|¬Gu|4α_↓|4v|¬G|4|¬W|4|πmax(|εu,|4v).|'
 2054  {A6}|πThese facts immediately suggest the following 
 2060  algorithm:|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m 
 2062  |≡B (|εBinary gcd algorithm).!|πGiven positive 
 2067  integers |εu |πand |εv, |πthis algorithm _nds 
 2074  their greatest common divisor.|'{A3}{I1.8H}|≡B|≡1|≡.|9[Find 
 2079  power of 2.] Set |εk|4|¬L|40, |πand then repeatedly 
 2087  set |εk|4|¬L|4k|4α+↓|41, u|4|¬L|4u/2, v|4|¬L|4v/2 
 2091  |πzero or more times until |εu |πand |εv |πare 
 2100  not both even.|'{A3}|≡B|≡2|≡.|9[Initialize.] 
 2104  (Now the original values of |εu |πand |εv |πhave 
 2113  been divided by |ε2|gk, |πand at least one of 
 2122  their present values is odd.) If |εu |πis odd, 
 2131  set |εt|4|¬L|4|→α_↓v |πand go to B4. Otherwise 
 2138  set |εt|4|¬L|4u.|'|π|≡B|≡3|≡.|9[Halve |εt.] (|πAt 
 2143  this point, |εt |πis even, and nonzero.) Set 
 2151  |εt|4|¬L|4t/2.|'{A3}|π|≡B|≡4|≡.|9[Is |εt |πeven?] 
 2155  If |εt |πis even, go back to B3.|'{A3}|≡B|≡5|≡.|9[Reset 
 2164  max(|εu,|4v).] |πIf |εt|4|¬T|40, |πset |εu|4|¬L|4t; 
 2169  |πotherwise set |εv|4|¬L|4|→α_↓t. |π(The larger 
 2174  of |εu |πand |εv |πhas been replaced by |ε|¬Gt|¬G, 
 2183  |πexcept perhaps during the _rst time this step 
 2191  is performed.)|'{A3}|≡B|≡6|≡.|9[Subtract.] Set 
 2195  |εt|4|¬L|4u|4α_↓|4v. |πIf |εt|4|=|↔6α=↓|40, |πgo 
 2199  back to B3. Otherwise the algorithm terminates 
 2206  with |εu|4|¬O|42|gk |πas the output.|'{A12}{IC}!|4|4|4|4As 
 2212  an example of Algorithm B, let us consider |εu|4α=↓|440902, 
 2221  v|4α=↓|424140, |πthe same numbers we have used 
 2228  to try out Euclid's algorithm. Step B1 sets |εk|4|¬L|41, 
 2237  u|4|¬L|420451, v|4|¬L|412070. |πThen |εt |πis 
 2242  set to |→α_↓12070, and replaced by |→α_↓6035; 
 2249  |εv |πis replaced by 6035, and the computation 
 2257  proceeds as follows:|'{A6}|∂!!!!!!|∂!!!!!!|∂!!!!!!!!!!!!!!!!
 2260  !!!|∂|E|;|>|εu|;v|;t|;>{A3}|>20451!|9|1|1|?6035!|4|4|4|4|?
 2269  !|→α+↓14416,|4|→α+↓7208,|4|→α+↓3604,|4|→α+↓1802,|4|→α+↓901;|
 2269  '>|>901!|9|1|1|?6035!|4|4|4|4|?!|→α_↓5134,|4|→α_↓2567;|'
 2275  >|>901!|9|1|1|?2567!|4|4|4|4|?!|→α_↓1666,|4|→α_↓833;|'
 2280  >|>901!|9|1|1|?833!|4|4|4|4|?!|→α+↓68,|4|→α+↓34,|4|→α+↓17;|'
 2285  >|>17!|9|1|1|?833!|4|4|4|4|?!|→α_↓816,|4|→α_↓408,|4|→α_↓204,
 2289  |4|→α_↓102,|4|→α_↓51;|'>|>17!|9|1|1|?51!|4|4|4|4|?
 2294  !|→α_↓34,|4|→α_↓17;|'>|>17!|9|1|1|?17!|4|4|4|4|?
 2299  !0.|'>{A6}|πThe answer is 17|4|¬O|42|g1|4α=↓|434. 
 2305  A few more iterations were necessary here than 
 2313  we needed with Algorithm A, but each iteration 
 2321  was somewhat simpler since no division steps 
 2328  were used.|'!|4|4|4|4A |¬m|¬i|¬x program for 
 2334  Algorithm B requires just a little more code 
 2342  than for Algorithm A. In order to make such a 
 2352  program fairly typical of a binary computer representation 
 2360  of Algorithm B, let us assume that |¬m|¬i|¬x 
 2368  is extended to include the following operators:|'
 2375  {A6}{J101}{H6}|*2⊃{H10}|9|¬s|¬l|¬b (shift left 
 2378  AX binary). C|4α=↓|46, F|4α=↓|46. The contents 
 2384  of registers A and X are ``shifted left'' M binary 
 2394  places; that is, |¬GrAX|¬G|4|¬L|4|¬G2|gMrAX|¬G 
 2398  mod |εB|g1|g0, |πwhere |εB |πis the byte size. 
 2406  (As with all |¬m|¬i|¬x shift commands, the signs 
 2414  of rA and rX are not a=ected.)|'{A3}{J101}{H6}|*2⊃{H10}|9|¬s|
 2421  ¬b|¬r (shift right AX binary). C|4α=↓|46, F|4α=↓|47. 
 2428  The contents of registers A and X are ``shifted 
 2437  right'' M binary places; that is, |¬GrAX|¬G|4|¬L|4|"l|¬GrAX|
 2443  ¬G/2|gM|"L.|'{A3}{J101}{H6}|*2⊃{H10}|9|¬j|¬a|¬e|¬, 
 2445  |¬j|¬a|¬o (jump A even, jump A odd). C|4α=↓|440; 
 2453  F|4α=↓|46, 7, respectively. A |¬j|¬m|¬p occurs 
 2459  if rA is even or odd, respectively.|'{A3}{J101}}h6}|*2⊃{H10}|
 2466  9|¬j|¬x|¬e|¬, |¬j|¬x|¬o (jump X even, jump X 
 2473  odd). C|4α=↓|447; F|4α=↓|46, 7, respectively. 
 2478  Analogous to |¬j|¬a|¬e|¬, |¬j|¬a|¬o.|'{A12}|≡P|≡r|≡o|≡g|≡r|≡
 2482  a|≡m |≡B (|εBineary gcd algorithm). |πAssume 
 2488  that |εu |πand |εv |πare single-precision, positive 
 2495  integers, stored respectively in locations |¬u 
 2501  and |¬v; this program uses Algorithm B to put 
 2510  gcd(|εu,|4v) |πinto rA. Register assignments: 
 2515  |εt|4|"o|4|πrA, |εk|4|"o|4|πrI1.|'{A6}|H!{U0}{H10L12M29}5832
folio 421 galley 4 WARNING: Much of this tape was unreadable!
 2517  0#Computer Programming!(A.-W./Knuth)!f. 421!ch 
 2520  4!g. 4|'{A12}{H9L11M29}|∂|*/|↔c|↔c!!|∂|\|π|¬a|¬b|¬s!!|∂|¬l|¬d
 2522  |¬a|¬n!!|∂|¬v|≤#|¬a|¬b|¬s|≤&!!|∂1|4α_↓|4B|4α+↓|4D!!|∂To|4B4|
 2522  4with|4|εt|4|¬L|4|→α_↓v|4|πif|4|εu|4|πis|4odd.|E|n|;
 2523  |L|ε|*/|↔c|↔O|L|\|π|¬a|¬b|¬s|L|¬e|¬q|¬u|L|¬1|¬.|¬5>
 2524  |L|ε|*/|↔c|↔P|\|π|L|¬b|¬1|L|¬e|¬n|¬t|¬1|L|¬0|L!!|1|1|11|L|εB|
 2524  */1|\.|9Find|4power|4of|4|*/2.>|L|↔c|↔L|\|π|L|L|¬l|¬d|¬x|L|¬u|
 2525  L!!|1|1|11|LrX|4|¬L|4|εu.>|L|*/|↔c|↔M|L|L|\|π|¬l|¬d|¬a|¬n|L|¬
 2526  v|L!!|1|1|11|LrA|4|¬L|4|→α_↓|εv.>|L|*/|↔c|↔C|L|L|\|π|¬j|¬m|¬p
 2527  |L|¬1|¬f|L!!|1|1|11>|L|ε|*/|↔c|↔o|\|π|L|¬2|¬h|L|¬s|¬r|¬b|L|¬1
 2528  |L|ε!!|1|1A|L|πHalve|4rA,|4rX.>|L|ε|*/|↔c|↔p|\|L|L|¬i|¬n|¬c|¬
 2529  1|L|¬1|L|ε!!|1|1A|Lk|4|¬L|4k|4α+↓|41.>|L|*/|↔c|↔l|\|L|L|¬s|¬t
 2530  |¬x|L|¬u|¬!!|1|1A|Lu|4|¬L|4u/2.>|L|*/|↔c|↔m|\|L|L|¬s|¬t|¬a|L|
 2531  ¬v|≤#|¬a|¬b|¬s|≤&|L!!|1|1A|Lv|4|¬L|4v/2.>*?{A12}|π|≡E|≡x|≡t|≡
 2532  e|≡n|≡s|≡i|≡o|≡n|≡s|≡.|9|4We can extend the methods 
 2537  used to calculate gcd(|εu,|4v) |πin order to 
 2544  solve some slightly more di∃cult problems. For 
 2551  example, assume that we want to compute the greatest 
 2560  common divisor of |εn |πintegers |εu|β1, u|β2,|4.|4.|4.|4,|4
 2566  u|βn.|'|π!|4|4|4One way to calculate gcd(|εu|β1, 
 2572  u|β2,|4.|4.|4.|4,|4u|βn), |πassuming that the 
 2576  |εu'|πs are all nonnegative, is to extend Euclid's 
 2584  algorithm in the following way: If all |εu|βj 
 2592  |πare zero, the greatest common divisor is taken 
 2600  to be zero; otherwise if only one |εu|βj |πis 
 2609  nonzero, it is the greatest common divisor; otherwise 
 2617  replace |εu|βk |πby |εu|βk |πmod |εu|βj |πfor 
 2624  all |εk|4|=|↔6α=↓|4j, |πwhere |εu|βj |πis the 
 2630  minimum of the nonzero |εu'|πs.|'!|4|4|4|4The 
 2636  algorithm sketched in the preceding paragraph 
 2642  is a natural generalization of Euclid's method, 
 2649  and it can be justi_ed in a similar manner. But 
 2659  there is a simpler method available, based on 
 2667  the easily veri_ed identity|'{A6}|πgcd(|εu|β1,|4u|β2,|4.|4.|
 2671  4.|4,|4u|βn)|4α=↓|4|πgcd{H12}({H10}|εu|β1,|4|πgcd(|εu|β2,|4.
 2671  |4.|4.|4,|4u|βn){H12}){H10}.|J!(14)|;{A6}|πTo 
 2673  calculate gcd(|εu|β1,|4u|β2,|4.|4.|4.|4,|4u|βn), 
 2675  |πwe may proceed as follows:|'{A6}|≡D|≡1|≡.|9Set 
 2681  |εd|4|¬L|4u|βn, j|4|¬L|4n|4α_↓|41.|'{A3}|π{I1.8H}|≡D|≡2|≡.|9
 2683  If |εd|4|=|↔6α=↓|41 |πand |εj|4|¬Q|40, |πset 
 2688  |εd|4|¬L|4|πgcd(|εu|βj,|4d) |πand |εj|4|¬L|4j|4α_↓|41 
 2691  |πand repeat this step. Otherwise |εd|4α=↓|4|πgcd(|εu|β1,|4.
 2696  |4.|4.|4,|4u|βn).|'{IC}{A12}|πThis method reduces 
 2700  the calculation of gcd|ε(u|β1,|4.|4.|4.|4,|4u|βn) 
 2704  |πto repeated calculations of the greatest common 
 2711  divisor of two numbers at a time. It makes use 
 2721  of the fact that gcd(|εu|β1,|4.|4.|4.|4,|4u|βj,|41)|4α=↓|41;
 2725   |πand this will be helpful, since we will already 
 2735  have gcd(|εu|βn|βα_↓|β1,|4u|βn)|4α=↓|41 |πover 
 2738  60 percent of the time if |εu|βn|βα_↓|β1 |πand 
 2746  |εu|βn |πare chosen at random. In most cases, 
 2754  the value of |εd |πwill rapidly decrease in the 
 2763  _rst few stages of the calculation, and this 
 2771  makes the remainder of the computation quite 
 2778  fast. Here Euclid's algorithm has an advantage 
 2785  over Algorithm B, in that its running time is 
 2794  primarily governed by the value of min(|εu,|4v), 
 2801  |πwhile the running time for Algorithm B is primarily 
 2810  governed by max(|εu,|4v); |πit would be reasonable 
 2817  to perform one iteration of Euclid's algorithm, 
 2824  replacing |εu |πby |εu |πmod |εv |πif |εu |πis 
 2833  much larger than |εv, |πand then to continue 
 2841  with Algorithm B.|'!|4|4|4|4The assertion that 
 2847  gcd(|εu|βn|βα_↓|β1,|4u|βn) |πwill be equal to 
 2852  unity more than 60 percent of the time for random 
 2862  inputs is a consequence of the following well-known 
 2870  result of number theory:|'{A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m 
 2875  |≡D (G. Lejeune Dirichlet, |εAbh. K|=4oniglich 
 2881  Preuss. Akad. Wiss. (1849), 69<83).|9|4If u and 
 2888  v are integers chosen at random, the probability 
 2896  that |πgcd(|εu,|4v)|4α=↓|41 is 6/|≤p|g2.|'|π!|4|4|4|4A 
 2901  precise formulation of this theorem, which carefully 
 2908  de_nes what is meant here by being ``chosen at 
 2917  random,'' appears in exercise 10 with a rigorous 
 2925  proof. Let us content ourselves here with a he{U0}{H10L12M29
folio 423 galley 5
 2933  }58320#Computer Programming!(A.-W./Knuth)!f. 
 2935  423!ch. 4!g. 5|'{A12}!|4|4|4|4If we assume, without 
 2942  proof, the existence of a well-de_ned probability 
 2949  |εp |πthat gcd(|εu,|4v) |πequals unity, then 
 2955  we can determine the probability that gcd(|εu,|4v)|4α=↓|4d 
 2962  |πfor any positive integer |εd; |πfor the latter 
 2970  event will happen only when |εu |πis a multiple 
 2979  of |εd, v |πis a multiple of |εd, |πand gcd(|εu/d,|4v/d)|4α=
 2988  ↓|41. |πThus the probability that gcd(|εu,|4v)|4α=↓|4d 
 2994  |πis equal to 1/|εd |πtimes 1/|εd |πtimes |εp, 
 3002  |πnamely |εp/d|g2. |πNow let us sum these probabilities 
 3010  over all possible values of |εd; |πwe should 
 3018  get|'{A6}|ε1|4α=↓|4|↔k|uc|)d|¬R1|)|4p/d|g2|4α=↓|4p(1|4α+↓|4|
 3019  f1|d34|)|4α+↓|4|f1|d39|)|4α+↓|4|f1|d316|)|4α+↓|1|1|¬O|4|¬O|4
 3019  |¬O).|;{A6}|πSince the sum 1|4α+↓|4|f1|d34|)|4α+↓|4|f1|d29|)
 3023  |1|1α+↓|4|¬O|4|¬O|4|¬O|4α=↓|4|εH|ur(2)|)|¬X|) 
 3024  |πis equal to |ε|≤p|g2 |πin order to make this 
 3033  equation come out right.!!|4|4|4|4Euclid's algorithm 
 3038  can be extended in another important way: We 
 3046  can calculate integers |εu|¬S |πand |εv|¬S |πsuch 
 3053  that|'{A6}|εuu|¬S|4α+↓|4vv|¬S|4α=↓|4|πgcd(|εu,|4v)|J!(15)|;
 3055  {A6}|πat the same time gcd(|εu,|4v) |πis being 
 3062  calculated. This extension of Euclid's algorithm 
 3068  can be described conveniently in vector notation:|'
 3075  {A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m |≡X (|εExtended 
 3078  Euclid's algorithm).|9|4|πGiven nonnegative integers 
 3082  |εu |πand |εv, |πthis algorithm determines a 
 3089  vector |ε(u|β1, u|β2, u|β3) |πsuch that |εuu|β1|4α+↓|4vu|β2|
 3095  4α=↓|4u|β3|4α=↓|4|πgcd(|εu,|4v). |πThe computation 
 3098  makes use of auxiliary vectors |ε(v|β1, v|β2, 
 3105  v|β3), (t|β1, t|β2, t|β3); |πall vectors are 
 3112  manipulated in such a way that the relations|'
 3120  {A6}|εut|β1|4α+↓|4vt|β2|4α=↓|4t|β3,!!uu|β1|4α+↓|4vu|β2|4α=↓|
 3120  4u|β3,!!uv|β1|4α+↓|4vv|β2|4α=↓|4v|β3|J!(16)|;
 3121  {A6}|πhold throughout the calculation.|'{A6}|≡X|≡1|≡.|9[Init
 3125  ialize.] Set |ε(u|β1, u|β2, u|β3)|4|¬L|4(1,|40,|4u), 
 3130  (v|β1,|4v|β2,|4v|β3)|4|¬L|4(0,|41,|4v).|'{A3}|π|≡X|≡2|≡.|9[I
 3131  s |εv|β3|4α=↓|40?] |πIf |εv|β3|4α=↓|40, |πthe 
 3136  algorithm terminates.|'{A3}|≡X|≡3|≡.|9[Divide, 
 3139  subtract.] Set |εq|4|¬L|4|"lu|β3/v|β3|"L, |πand 
 3143  then set|'{A6}|ε(t|β1,|4t|β2,|4t|β3)|4|¬L|4(u|β1,|4u|β2,|4u|
 3145  β3)|4α_↓|4(v|β1,|4v|β2,|4v|β3)q,|;{A6}(u|β1,|4u|β2,|4u|β3)|4
 3146  |¬L|4(v|β1,|4v|β2,|4v|β3),!!(v|β1,|4v|β2,|4v|β3)|4|¬L|4(t|β1
 3146  ,|4t|β2,|4t|β3).|;{A6}!|4|4|4|4|πReturn to step 
 3150  X2.|'{A12}!|4|4|4|4For example, let |εu|4α=↓|440902, 
 3155  v|4α=↓|424140. |πAt step X2 we have|'{A6}|∂!!!!|∂!!!!|∂!!!!|
 3161  ∂!!!!|∂!!!!|∂!!!!|∂!!!!|∂|E|;|>|εq|;u|β1|;u|β2|;
 3166  u|β3|;v|β1|;v|β2|;v|β3|;>|>#|;1!|1|?0|4|4|4|?
 3175  40902|9|1|1|?0|4|4|4|?1|4|4|4|?24140|9|1|1|?>
 3180  |>1|;0!|1|?1|4|4|4|4|?24140|9|1|1|?1|4|4|4|?|→α_↓1|4|4|4|?
 3187  16762|9|1|1|?>|>1|;1!|1|?|→α_↓1|4|4|4|?16762|9|1|1|?
 3194  |→α_↓1|4|4|4|?2|4|4|4|?7378|9|1|1|?>|>2|;|→α_↓1!|1|?
 3201  2|4|4|4|?7378|9|1|1|?3|4|4|4|?|→α_↓5|4|4|4|?2006|9|1|1|?
 3206  >|>3|;3!|1|?|→α_↓5|4|4|4|?2006|9|1|1|?|→α_↓10|4|4|4|117|4|4|
 3212  4|?1360|9|1|1|?>|>1|;|→α_↓10!|1|?17|4|4|4|?1360|9|1|1|?
 3220  13|4|4|4|?|→α_↓22|4|4|4|?646|9|1|1|?>|>2|;13!|1|?
 3227  |→α_↓22|4|4|4|?646|9|1|1|?|→α_↓36|4|4|4|?61|4|4|4|?
 3231  68|9|1|1|?>|>9|;|→α_↓36!|1|?61|4|4|4|?68|9|1|1|?
 3238  337|4|4|4|?|→α_↓571|4|4|4|?34|9|1|1|?>|>2|;337!|1|?
 3245  |→α_↓571|4|4|4|?34|9|1|1|?|→α_↓710|4|4|4|?1203|4|4|4|?
 3249  0|9|1|1|?>{A6}|πTherefore the solution is 337|4|¬O|440902|4α
 3255  _↓|4571|4|¬O|424140|4α=↓|434|4α=↓|4gcd(40902,|424140).|'
 3256  {A6}!|4|4|4|4The validity of Algorithm X follows 
 3262  from (16) and the fact that the algorithm is 
 3271  identical to Algorithm A with respect to its 
 3279  manipulation of |εu|β3 |πand |εv|β3. |πA detailed 
 3286  proof of Algorithm X is discussed in Section 
 3294  1.2.1. Gordon H. Bradley has observed that we 
 3302  can avoid a good deal of the calculation in Algorithm 
 3312  X by suppressing |εu|β1, v|β1, |πand |εt|β1; 
 3319  |πthen |εu|β1 |πcan be determined afterwards 
 3325  using the relation |εuu|β1|4α+↓|4vu|β2|4α=↓|4u|β3.|'
 3329  !|4|4|4Exercise 14 shows that the values of |¬G|εu|β1|¬G, 
 3337  |¬Gu|β2|¬G, |¬Gv|β1|¬G, |¬Gv|β2|¬G |πremain bounded 
 3342  by the size of the inputs |εu |πand |εv.|'|π!|4|4|4|4Algorit
 3351  hm B, which computes the greatest common divisor 
 3359  properties of binary notation, does |εnot |πextend 
 3366  in a similar way to an algorithm that solves 
 3375  (15). For some instructive extensions to Algorithm 
 3382  X, see exercises 18 and 19 in Section 4.6.1.|'
 3391  {A6}!|4|4|4|4The ideas underlying Euclid's algorithm 
 3396  can also be applied to _nd a |εgeneral solution 
 3405  in integers |πof any set of linear equations 
 3413  with integer coe∃cients. For example, suppose 
 3419  that we want to _nd all integers |εw, x, y, z 
 3430  |πwhich satisfy the two equations|'{A6}|ε|∂10w|4α+↓|43x|4α+↓
 3435  |43y|4|∂α+↓|48z|4α=↓|41,|J!(17)|;{A4}|L|9|16w|4α_↓|47x|Lα_↓|
 3436  45z|4α=↓|42.|J!(18)>{A6}|πWe can introduce a 
 3441  new variable|'{A6}|"l10/3|"L|εw|4α+↓|4|"l3/3|"Lx|4α+↓|4|"l3/
 3443  3|"Ly|4α+↓|4|"l8/3|"Lz|4α=↓|43w|4α+↓|4x|4α+↓|4y|4α+↓|42z|4α=
 3443  ↓|4t|β1,|;{A6}|πand use it to eliminate |εy; 
 3450  |πEq. (17) becomes|'{A6}(10|4mod|43)|εw|4α+↓|4(3|4|πmod|43)|
 3453  εx|4α+↓|43t|β1|4α+↓|4(8|4|πmod|43)|εz|4α=↓|4w|4α+↓|43t|β1|4α
 3453  +↓|42z|4α=↓|41,|J!(19)|;{A6}|πand Eq. (18) remains 
 3458  unchanged. The new equation (19) may be used 
 3466  to eliminate |εw, |πand (18) becomes|'{A6}|ε6(1|4α_↓|43t|β1|
 3472  4α_↓|42z)|4α_↓|47x|4α_↓|45z|4α=↓|42;|;{A6}|πthat 
 3474  is,|'{A6}|ε7x|4α+↓|418t|β1|4α+↓|417z|4α=↓|44.|J!(20)|;
 3476  {A6}|πNow as before we introduce a new variable|'
 3484  {A6}|εx|4α+↓|42t|β1|4α+↓|42z|4α=↓|4t|β2|;{A6}|πand 
 3486  eliminate |εx |πfrom (20):|'{A6}|ε7t|β2|4α+↓|44t|β1|4α+↓|43z
 3490  |4α=↓|44.|J!(21)|;{A6}|πAnother new variable 
 3494  can be introduced in the same fashion, in order 
 3503  to eliminate the variable |εz, |πwhich has the 
 3511  smallest coe∃cient:|'{A6}|ε2t|β2|4α+↓|4t|β1|4α+↓|4z|4α=↓|4t|
 3513  β3.|;{A6}|πEliminating |εz |πfrom (21) yields|'
 3519  {A6}|εt|β2|4α+↓|4t|β1|4α+↓|43t|β3|4α=↓|44,|J!(22)|;
 3520  {A6}|πand this equation, _nally, can be used 
 3527  to eliminate |εt|β2. |πWe are left with two independent 
 3536  variables, |εt|β1 |πand |εt|β3; |πsubstituting 
 3541  back for the original variables, we obtain the 
 3549  general solution|'{A6}|εw|4|∂α=↓|4|4|4|4|417|4α_↓|4|9|15t|β1
 3551  |4α_↓|414t|β3,|;{A4}| x|4|Lα=↓|4|4|4|4|420|4α_↓|4|9|15t|β1|4
 3552  α_↓|417t|β3,>{A4}| y|4|Lα=↓|4|→α_↓55|4α+↓|419t|β1|4α+↓|445t|
 3553  β3,>{A4}| z|4|Lα=↓|4|9|1|→α_↓8|4α+↓|4|4|4|4|4|4t|β1|4α+↓|4|9
 3554  |17t|β3.>{B24}(23)|?{A24}{A6}|πIn other words, 
 3559  all integer solutions |ε(w,|4x,|4y,|4z) |πto 
 3564  the original equations (17), (18) are obtained 
 3571  from (23) by letting |εt|β1 |πand |εt|β3 |πindependently 
 3579  run through all integers.|'!|4|4|4|4The general 
 3585  method which has just been illustrated is based 
 3593  on the following procedure: Find a nonzero coe∃cient 
 3601  |εc |πof smallest absolute value in the system 
 3609  of equations. Suppose that this c{U0}{H10L12M29}59320#Comput
folio 426 galley 6
 3614  er Programming!(A.-W./Knuth)!ch. 4!f. 426!g. 
 3618  6|'{A12}|πwe may assume for simplicity that |εc|4|¬Q|40. 
 3626  |πIf |εc|4α=↓|41, |πuse this equation to eliminate 
 3633  the variable |εx|β0 |πfrom the other equations 
 3640  remaining in the system; then repeat the procedure 
 3648  on the remaining equations. (If no more equations 
 3656  remain, the computation stops, and a general 
 3663  solution in terms of the variables not yet eliminated 
 3672  has essentially been obtained.) If |εc|4|¬Q|41, 
 3678  |πthen if |εc|β1 |πmod |εc|4α=↓|1|1|¬O|4|¬O|4|¬O|1|1α=↓|4c|β
 3682  k |πmod |εc|4α=↓|40 |πwe must have |εd |πmod 
 3690  |εc|4α=↓|40, |πotherwise there is no integer 
 3696  solution; divide both sides of (24) by |εc |πand 
 3705  eliminate |εx|β0 |πas in the case |εc|4α=↓|41. 
 3712  |πFinally, if |εc|4|¬Q|41 |πand not all of |εc|β1 
 3720  |πmod |εc,|4.|4.|4.|4,|4c|βk |πmod |εc |πare 
 3725  zero, then introduce a new variable.|'{A6}|ε|"lc/c|"Lx|β0|4α
 3731  +↓|4|"lc|β1/c|"Lx|β1|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|"lc|βk/c
 3731  |"Lx|βk|4α=↓|4t;|J!(25)>{A6}|πeliminate the variable 
 3735  |εx|β0 |πfrom the other equations, in favor of 
 3743  |εt, |πand replace the original equation (24) 
 3750  by|'{A6}|εct|4α+↓|4(c|β1|4|πmod|4|εc)x|β1|4α+↓|1|1|¬O|4|¬O|4
 3751  |¬O|1|1α+↓|4(c|βk|4|πmod|4|εc)x|βk|4α=↓|4d.|J!(26)|;
 3752  {A6}|π(Cf. (19) and (21) in the above example.)|'
 3760  !|4|4|4|4This process must terminate, since each 
 3766  step either reduces the number of equations or 
 3774  the size of the smallest nonzero coe∃cient in 
 3782  the system. A study of the above procedure will 
 3791  reveal its intimate connection with Euclid's 
 3797  algorithm. The method is a comparatively simple 
 3804  means of solving linear equations when the variables 
 3812  are required to take on integer values only.|'
 3820  {A12}|≡H|≡i|≡g|≡h|≡-|≡p|≡r|≡e|≡c|≡i|≡s|≡i|≡o|≡n 
 3821  |≡c|≡a|≡l|≡c|≡u|≡l|≡a|≡t|≡i|≡o|≡n|≡.|9|4If |εu 
 3823  |πand |εv |πare very large integers, requiring 
 3830  a multiple-precision representation, the binary 
 3835  method (Algorithm B) is a simple and reasonably 
 3843  e∃cient means of calculating their greatest common 
 3850  divisor.|'!|4|4|4|4By contrast, Euclid's algorithm 
 3855  seems much less attractive, since step A2 requires 
 3863  a multiple precision division of |εu |πby |εv. 
 3871  |πBut this di∃culty is not really as bad as it 
 3881  seems, since we will prove in Section 4.5.3 that 
 3890  the quotient |"l|εu/v|"L |πis almost always very 
 3897  small; for example, assuming random inputs, the 
 3904  quotient |ε|"lu/v|"L |πwill be less than 1000 
 3911  approximately 99.856 percent of the time. Therefore 
 3918  it is almost always possible to _nd |ε|"lu/v|"L 
 3926  |πand (|εu |πmod |εv) |πusing single precision 
 3933  calculations, together with the comparatively 
 3938  simple operation of calculating |εu|4α_↓|4qv 
 3943  |πwhere |εq |πis a single-precision number.|'
 3949  !|4|4|4|4A signi_cant improvement in the speed 
 3955  of Euclid's algorithm when high-precision numbers 
 3961  are involved can be achieved by using a methmo*?*?-precision 
 3970  numbers are involved can be achieved by using 
 3978  a method due to D. H. Lehmar [|εAMM |≡4|≡5 (1938), 
 3988  227<233]. |πWorking only with the leading digits 
 3995  of large numbers, it is possible to do most of 
 4005  the calculations with single-precision arithmetic, 
 4010  and to make a substantial reduction in the number 
 4019  of multiple-precision operations involved.|'!|4|4|4|4Lehmer'
 4023  s method can be illustrated on the eight-digit 
 4031  numbers |εu|4α=↓|427182818, v|4α=↓|410000000, 
 4034  |πassuming that we are using a machine with only 
 4043  four-digit words. Let |εu|¬S|4α=↓|42718, v|¬S|4α=↓|41001, 
 4048  u|¬C|4α=↓|42719, v|¬C|4α=↓|41000; |πthen |εu|¬S/v|¬S 
 4052  |πand |εu|¬C/v|¬C |πare approximations to |εu/v, 
 4058  |πwith|'{A6}|εu|¬S/v|¬S|4|¬W|4u/v|4|¬W|4u|¬C/v|¬C.|J!(27)|;
 4060  {A6}|πThe ratio |εu/v |πdetermines the sequence 
 4066  of quotients obtained in Euclid's algorithm; 
 4072  if we carry out Euclid's algorithm simultaneously 
 4079  on the single-precision values (|εu|¬S,|4v|¬S) 
 4084  |πand |ε(u|¬C,|4v|¬C) |πuntil we get a di=erent 
 4091  quotient, it is not di∃cult to see that the same 
 4101  sequence of quotients would have appeared to 
 4108  this point if we had worked with the multiple 
 4117  precision numbers |ε(u,|4v). |πThus consider 
 4122  what happens when Euclid's algorithm is applied 
 4129  to |ε(u|¬S,|4v|¬S) |πand to |ε(u|¬C,|4v|¬C):|'
 4134  {A6}|∂!!!!|∂!!!!|∂!!!!|∂!!!!|∂!!!!|∂!!!!|∂!!!!|∂|E|;
 4135  |>u|¬S|;v|¬S|;q|¬S|;|;u|¬C|;v|¬C|;q|¬C|;>{A3}|>
 4145  2718|;1001|;|9|12|;|;2719|;1000|;|9|12|;>|>1001|;
 4155  |9|1716|;|9|11|;|;10000|;|9|1719|;1|;>|>|9|1716|;
 4164  |9|1285|;|9|12|;|;|9|1719|;|9|1281|;2|;>|>|9|1285|;
 4173  |9|1146|;|9|11|;|;|9|1281|;|9|1157|;1|;>|>|9|1146|;
 4182  |9|1139|;|9|11|;|;|9|1157|;|9|1124|;3|;>|>|9|1139|;
 4191  !|4|4|47|;19|;|;|9|1124|;!|1|133|;3|;>{A6}|πAfter 
 4199  six steps we _nd that |εq|¬S|4|=|↔6α=↓|4q|¬C, 
 4205  |πso the single-precision calculations are suspended; 
 4211  we have gained the knowledge that the calculation 
 4219  would have proceeded as follows if we had been 
 4228  working with the original multiple-precision 
 4233  numbers:|'{A6}|h|ε|→α_↓4u|β0|4|∂α+↓|411v|β0!!!!|→α_↓4u|β0|4|
 4234  ∂α+↓|411v|β0!!!!1|E|n|;|Lu|Lv|Lq>|Lu|β0|Lv|β0|L2>
 4237  |Lv|β0| u|β0|4|Lα_↓|42v|β0|L1>| u|β0|4|Lα_↓|42v|β0| |→α_↓u|β
 4238  0|4|Lα+↓|43v|β0|L2|J!(28)>| |→α_↓u|β0|4|Lα+↓|43v|β0| 3u|β0|4
 4239  |Lα_↓|48v|β0|L1>| 3u|β0|4|Lα_↓|L8v|β0| |→α_↓4u|β0|4|Lα+↓|L11
 4240  v|β0|L1>| |→α_↓4u|β0|4|Lα+↓|411v|β0| 7u|β0|4|Lα_↓|419v|β0|L?
 4241  >{A6}|π(The next quotient lies somewhere between 
 4248  3 and 19.) No matter how many digits are in |εu 
 4259  |πand |εv, |πthe _rst _ve steps of Euclid's algorithm 
 4268  would be the same as (28), so long as (27) holds. 
 4279  We can there fore avoid the multiple-precision 
 4286  operations of the _rst _ve steps, and replace 
 4294  them all by a multiple-precision calculation 
 4300  of |→α_↓4|εu|β0|4α+↓|411v|β0 |πand |ε7u|β0|4α_↓|419v|β0. 
 4304  |πIn this case we would obtain |εu|4α=↓|41268728, 
 4311  v|4α=↓|4279726; |πthe calculation could now proceed 
 4317  with |εu|¬S|4α=↓|41268, v|¬S|4α=↓|4280, u|¬C|4α=↓|41269, 
 4321  v|¬C|4α=↓|4279, |πetc. With a larger accumulator, 
 4327  more steps could be done by single-precision 
 4334  calculations; our example showed that only _ve 
 4341  cycles of Euclid's algorithm were combined into 
 4348  one multiple step, but with (say) a word size 
 4357  of 10 digits we could do about twelve cycles 
 4366  at a time. (Results proved in Section 4.5.3 imply 
 4375  that the number of multiple-precision cycles 
 4381  which can be replaced at each iteration is essentially 
 4390  proportional to the logarithm of the word size 
 4398  used in the single-precision calculations.)|'
 4403  !|4|4|4|4Lehmer's method can be formulated as 
 4409  follows:|'{A12}|≡A|≡l|≡g|≡o|≡r|≡i|≡t|≡h|≡m |≡L 
 4412  (|εEuclid's algorithm for large numbers)|9|4|πLet 
 4417  |εu, v |πbe nonnegative integers, with |εu|4|¬R|4v, 
 4424  |πrepresented in multiple precision. This algorithm 
 4430  computes the greatest common divisor of |εu |πand 
 4438  |εv, |πmaking use of auxiliary single-precision 
 4444  |εp-|πdigit variables |εu, v, A, B, C, D, T, 
 4453  |πand |εq, |πand auxiliary multiprecision variables 
 4459  |εt |πand |εw.|'{A12}|π{I1.8H}|≡L|≡1|≡.|9[Initialize.] 
 4463  If |εv |πis small enough to be represented as 
 4472  a single-precision value, calculate gcd(|εu,|4v) 
 4477  |πby Algorithm A and terminate the computation. 
 4484  Otherwise, let |ε|=7u |πbe the |εp |πleading 
 4491  digits of |εu, |πand let |ε|=7v |πbe the |εp 
 4500  |πcorresponding digits of |εv; |πin other words, 
 4507  if radix |εb |πnotation is being used, |ε|=7u|1|1|¬L|1|1|¬lu
 4514  /b|gk|"L, |=#v|4|¬L|4|"lv/b|gk|"L, |πwhere |εk 
 4518  |πis as small as possible consistent with the 
 4526  condition |ε|=7u|4|¬W|4b|gp.|'!!!|4|4|4|4|πSet 
 4529  |εA|4|¬L|41, B|4|¬L|40, C|4|¬L|40, D|4|¬L|41. 
 4533  (|πThese variables represent the coe∃cients in 
 4539  (28), where|'{A6}!!|εu|4α=↓|4Au|β0|4α+↓|4Bv|β0,!!v|4α=↓|4Cu|
 4541  β0|4α+↓|4Dv|β0|J!(29)|;{A6}|π!!in the equivalent 
 4545  actions of Algorithms A on multiprecision numbers. 
 4552  We also have|'{A6}!!|εu|¬S|4α=↓|4|=7u|4α+↓|4B,!!v|¬S|4α=↓|4|
 4555  =7v|4α+↓|4D,!!u|¬C|4α=↓|4|=7u|4α+↓|4A,!!v|¬C|4α=↓|4|=7v|4α+↓
 4555  |4C|J!(30)|;{A6}|π!!in terms of the notation 
 4561  in the example worked above.)|'{A3}|≡L|≡2|≡.|9[Test 
 4567  quotient.] Set |εq|4|¬L|4|"l(|=7u|4α+↓|4A)/(|=7v|4α+↓|4C)|"L
 4569  . |πIf |εq|4|=|↔6α=↓|4|"l(|=7u|4α+↓|4B)/(|=7v|4α+↓|4D)|"L, 
 4572  |πgo to step L4. (This step tests if |εq|¬S|4|=|↔6α=↓|4q|¬C 
 4581  |πin the notation of the above example. Note 
 4589  that single-precision over⊗ow can occur in special 
 4596  circumstances during the computation in this 
 4602  step, but only when |ε|=7u|4α=↓|4b|gp|4α_↓|41 
 4607  |πand |εA|4α=↓|41 |πor when |ε|=7v|4α=↓|4b|gp|4α_↓|41 
 4612  |πand |εD|4α=↓|41; |πthe conditions|'|'a6}|ε!!0|4|¬E|4|=7u|4
 4617  α+↓|4A|4|¬E|4b|gp,!!0|4|¬E|4|=7v|4α+↓|4{U0}{H10L12M29}58320#
folio 429 galley 7
 4617  Computer Programming!(A.-W./Knuth)!ch. 4f. 429!g. 
 4621  7|'{A12}{I1.8H}will always hold, because of (30). 
 4628  It is possible to have |ε|=7v|4α+↓|4C|4α=↓|40 
 4634  |πor |ε|=7v|4α+↓|4D|4α=↓|40, |πbut not both simultaneously; 
 4640  therefore division by zero in this step is taken 
 4649  to mean ``Go directly to L4.'')|'{A3}|π|≡L|≡3|≡.|9[Emulate 
 4656  Euclid.] Set |εT|4|¬L|4A|4α_↓|4qC, A|4|¬L|4C, 
 4660  C|4|¬L|4T, T|4|¬L|4B|4α_↓|4qD, B|4|¬L|4D, D|4|¬L|4T, 
 4664  T|4|¬L|4|=7u|4α_↓|4q|=7v, |=7u|4|¬L|4|=7v, |=7v|4|¬L|4T, 
 4667  |πand go back to step L2. (These single-precision 
 4675  calculations are the equivalent of multiple-precision 
 4681  operations, as in (28), under the conventions 
 4688  of (29).)|'{A3}|≡L|≡4|≡.|9[Multiprecision step.] 
 4692  If |εB|4α=↓|40, |πset |εt|4|¬L|4u |πmod |εv, 
 4698  u|4|¬L|4v, v|4|¬L|4t |πusing multiple-precision 
 4702  division. (This happens only if the single-precision 
 4709  operations cannot simulate any of the multiprecision 
 4716  ones. It implies that Euclid's algorithm requires 
 4723  a very large quotient, and this is an extremely 
 4732  rare occurrence.) Otherwise, set |εt|4|¬L|4Au, 
 4737  t|4|¬L|4t|4α+↓|4Bv, w|4|¬L|4Cu, w|4|¬L|4w|4α+↓|4Dv, 
 4740  u|4|¬L|4t, v|4|¬L|4w (|πusing straightforward 
 4744  multiprecision operations). Go back to step L1.|'
 4751  {A6}{IC}!|4|4|4|4The values of |εA, B, C, D |πremain 
 4759  as single-precision numbers throughout this calculation, 
 4765  because of (31).|'!|4|4|4|4Algorithm L requires 
 4771  a somewhat complicated program than Algorithm 
 4777  B, but with large numbers it will be faster on 
 4787  many computers. Algorithm B can, however, be 
 4794  speeded up in a similar way (see exercise 34), 
 4803  to the point where it continues to win. Algorithm 
 4812  L has the advantage that it can be extended, 
 4821  as in Algorithm X (see exercise 17); furthermore, 
 4829  it determines the sequence of quotients obtained 
 4836  in Euclid's algorithm, and this yields the regular 
 4844  continued fraction expansion of a real number 
 4851  (see exercise 4.5.3<18).|'{A12}|≡A|≡n|≡a|≡l|≡y|≡s|≡i|≡s 
 4855  |≡o|≡f |≡t|≡h|≡e |≡b|≡i|≡n|≡a|≡r|≡y |≡a|≡l|≡g|≡o|≡r|≡i|≡t|≡h
 4858  |≡m|≡.|9|4Let us conclude this section by studying 
 4865  the running time of Algorithm B, in order to 
 4874  justify the formulas stated earlier.|'!|4|4|4|4An 
 4880  exact determination of the behavior of Algorithm 
 4887  B appears to be exceedingly di∃cult to derive, 
 4895  but we can begin to study it by means of an approximate 
 4907  model of its behavior. Suppose that |εu |πand 
 4915  |εv |πare odd numbers, with |εu|4|¬Q|4v |πand|'
 4922  {A6}|ε|"l|πlg|β2|4|εu|"L|4α=↓|4m,!!|π|"llg|β2|4|εv|"L|4α=↓|4
 4922  n.|Ja(32)|;{A6}|π{H12}({H10}Thus, |εu |πis an 
 4927  (|εm|4α+↓|41)-|πbit number, and |εv |πis an (|εn|4α+↓|41)-|π
 4933  bit number.{H12}){H10}. Algorithm B forms |εu|4α_↓|4v 
 4939  |πand shifts this quantity right until obtaining 
 4946  an odd number |εu|¬S |πwhich replaces |εu. |πUnder 
 4954  random conditions, we would expect to have|'{A6}|εu|¬S|4α=↓|
 4961  4(u|4α_↓|4v)/2|;{A6}|πabout one-half of the time,|'
 4967  {A6}|εu|¬S|4α=↓|4(u|4α_↓|4v)/4|;{A6}|πabout one-fourth 
 4970  of the time,|'{A6}|εu|¬S|4α=↓|4(u|4α_↓|4v)/8|;
 4974  {A6}|πabout one-eighth of the time, and so on. 
 4982  We have|'{A6}|ε|"l|πlg|β2|4|εu|¬S|"L|4α=↓|4m|4α_↓|4k|4α_↓|4r
 4984  ,|J!(33)|;{A6}|πwhere |εk |πis the number of 
 4991  places that |εu|4α_↓|4v |πis shifted right, and 
 4998  where |εr |πis |"llg|β2|4|εu|"L|4α_↓|4|"l|πlg|β2|4(|εu|4α_↓|
 5001  4v)|"L, |πthe number of bits lost at the left 
 5010  during the subtraction of |εv |πfrom |εu. |πNote 
 5018  that |εr|4|¬E|41 |πwhen |εm|4|¬R|4n|4α+↓|42, 
 5022  |πand |εr|4|¬R|41 |πwhen |εm|4α=↓|4n. |πFor simplicity, 
 5028  we will assume that |εr|4α=↓|40 |πwhen |εm|4|=|↔6α=↓|4n 
 5035  |πand that |εr|4α=↓|41 |πwhen |εm|4α=↓|4n, |πalthough 
 5041  this lower bound tends to make |εu|¬S |πseem 
 5049  larger than it usually is.|'!|4|4|4|4The approximate 
 5056  model we shall use to study Algorithm B is based 
 5066  solely on the values |εm|4α=↓|4|"l|πlg|β2|4|εu|"L 
 5071  |πand |εn|4α=↓|4|"l|πlg|β2|4|εv|"L |πthroughout 
 5074  the course of the algorithm, not on the actual 
 5083  values of |εu |πand |εv. |πLet us call this approximation 
 5093  a |εlattice-point model, |πsince we will say 
 5100  that we are ``at the point (|εm, n)'' |πwhen 
 5109  |"llg|β2|4|εu|"L|4α=↓|4m |πand |"l|πlg|β2|4|εv|"L|4α=↓|4n. 
 5112  |πFrom point |ε(m,|4n) |πthe algorithm takes 
 5118  us to |ε(m|¬S,|4n) |πif |εu|4|¬Q|4v, |πor to 
 5125  |ε(m,|4n|¬S) |πif |εu|4|¬W|4v, |πor terminates 
 5130  if |εu|4α=↓|4v. |πFor example, the calculation 
 5136  starting with |εu|4α=↓|420451, v|4α=↓|46035 |πwhich 
 5141  is tabulated after Algorithm B begins at the 
 5149  point (14, 12), then goes to (9, 12), (9,|411), 
 5158  (9,|49), (4,|49), (4,|45), (4,|44), and terminates. 
 5164  In line with the comments of the preceding paragraph, 
 5173  we will make the following assumptions about 
 5180  the probability that we reach a given point just 
 5189  after point |ε(m,|4n):|'{A6}|πCase 1, |εm|4|¬Q|4n.|;
 5195  {A3}|∂!!!!!!|∂!!!!!!|∂|E|;|π|>Next|4point|;Probability|;
 5199  >{A3}|>|ε(m|4α_↓|41,|4n)|;|f1|d32|)|;>|>(m|4α_↓|42,|4n)|;
 5206  |f1|d34|)|;>|>|¬O|4|¬O|4|¬O|;|¬O|4|¬O|4|¬O|;>
 5212  |>(1,|4n)|;|f1|d32|)|1|gm|gα_↓|g1|;>|>(0,|4n)|;
 5218  |f1|d32|)|1|gm|gα_↓|g1|;>{A12}|πCase 2, |εm|4|¬W|4n.|;
 5223  {A3}|π|>Next|4point|;Probability|;>{A3}|>|ε(m,|4n|4α_↓|41)|;
 5229  |f1|d32|)|;>|>(m,|4n|4α_↓|42)|;|f1|d34|)|;>|>
 5236  |¬O|4|¬O|4|¬O|;|¬O|4|¬O|4|¬O|;>|>(m,|41)|;|f1|d32|)|1|gn|gα_
 5241  ↓|g1|;>|>(m,|40)|;|f1|d32|)|1|gn|gα_↓|g1|;>{A12}|πCase 
 5248  3, |εm|4α=↓|4n|4|¬Q|40.|;{A3}|∂!!!!!!!!!|∂!!!!!!|∂|E|;
 5251  |>|πNext|4point|;Probability|;>{A3}|>|ε(m|4α_↓|42,|4n),|4(m,
 5256  |4n|4α_↓|42)|;|f1|d34|),|4|f1|d34|)|;>|>(m|4α_↓|43,|4n),|4(m
 5260  ,|4n|4α_↓|43)|;|f1|d38|),|4|f1|d38|)|;>|>|¬O|4|¬O|4|¬O|;
 5265  |¬O|4|¬O|4|¬O|;>|>(0,|4n),|4(m,|40)|;|f1|d32|)|1|gm,|4|f1|d3
 5269  2|)|1|gm|;>|>|πterminate|;|ε|f1|d32|)|1|gm|gα_↓|g1|;
 5274  >{A12}|πFor example, from points (5,|43) the 
 5281  lattice-point model would go to points (4,|43), 
 5288  (3,|43), (2,|43), (1,|43), (0,|43) with the respective 
 5295  probabilities |f1|d32|), |f1|d34|), |f1|d38|), 
 5299  |f1|d316|), |f1|d316|); from (4,|44) it would 
 5305  go to (2,|44), (1,|44), (0,|44), (4,|42), (4,|41), 
 5312  (4,|40), or would terminate, with the respective 
 5319  probabilities |f1|d34|), |f1|d38|), |f1|d316|), 
 5323  |f1|d34|), 1|d38|), |f1|d316|), |f1|d38|). When 
 5328  |εm|4α=↓|4n|4α=↓|40, |πthe formulas above do 
 5333  not apply; the algorithm always terminates in 
 5340  such a case, since |εm|4α=↓|4n|4α=↓|40 |πimplies 
 5346  that |εu|4α=↓|4v|4α=↓|41.|'!|4|4|4|4|πThe detailed 
 5350  calculations of exercise 18 show that this lattice-point 
 5358  model is somewhat pessimistic. In fact, when 
 5365  |εm|4|¬Q|43 |πthe actual probability that Algorithm 
 5371  B goes from |ε(m,|4m) |πto one of the two points 
 5381  (|εm|4α_↓|42, m) |πor |ε(m, m|4α_↓|42) |πis equal 
 5388  to |f1|d38|), although we have assumed the value 
 5396  |f1|d32|); the algorithm actually goes from |ε(m,|4m) 
 5403  |πto |ε(m|4α_↓|43,|4m) |πor |ε(m,|4m|4α_↓|43) 
 5407  |πwith probability |f7|d332|), not |f1|d34|); 
 5412  it actually goes from |ε(m|4α+↓|41, m) |πto |ε(m,|4m) 
 5420  |πwith probability |f1|d38|), not |f1|d32. The 
 5426  probabilities in the model are nearly correct 
 5433  when |ε|¬Gm|4α_↓|4n|¬G |πis large, but when |ε|¬Gm|4α_↓|4n|¬
 5439  G |πis small the model predicts slower convergence 
 5447  than is actually obtained. In spite of the fact 
 5456  that our model is not a completely faithful representation 
 5465  of Algorithm B, it has one great virtue, namely 
 5474  that it can be completely analyzed*3 Furthermore, 
 5481  empirical experiments with Algorithm B show that 
 5488  the behavior predicted by the lattice-point model 
 5495  is analogous to the true behavior.|'{U0}{H10L12M29}|π!|4|4|4
 5501  |4An analysis of the lattice-point model can 
 5508  be carried out by solving the following rather 
 5516  complicated set of double recurrence relations:|'
 5522  |h|εA|βm|βn|4α=↓|4c|4α+↓|4|9A|β(|βm|βα_↓|β1|β)|βn|4α+↓|1|1|¬
 5522  O|4|≡O|4|≡O|1|1α+↓|42|gm|gα_↓|g1|4A|β1|βn|4α+↓|42|gm|gα_↓|g1
 5522  |4A|β0|βn,!|∂|πif!|∂just|E|n|'{A6}|εA|βm|βm|4α=↓|4a|4α+↓|4|f
 5523  1|d32|)A|βm|β(|βm|βα_↓|β2|β)|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|
 5523  (1|d22|gm|gα_↓|g1|)|4A|βm|β0|4α+↓|4|(b|d22|gm|gα_↓|g1|)|4,|L
 5523  |πif|L|εm|4|¬R|41;>{A6}A|βm|βn|4α=↓|4c|4α+↓|4|f1|d32|)A|β(|β
 5524  m|βα_↓|β1|β)|βn|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|(1|d22|gm|gα_
 5524  ↓|g1|)|4A|β1|βn|4α+↓|4|(1|d22|gm|gα_↓|g1|)|4A|β0|βn,|L|πif|L
 5524  |εm|4|¬Q|4n|4|¬R|40;>{A6}A|βm|βn|4α=↓|4A|βn|βm,|L|πif|L|εn|4
 5525  |¬Q|4m|4|¬R|40.>{a3}(34)|?{A6}|πThe problem is 
 5530  to solve for |εA|βm|βn |πin terms of |εm, n, 
 5539  |πand the parameters |εa, b, c |πand |εA|β0|β0. 
 5547  |πThis is an interesting set of recurrence equations, 
 5555  which have an interesting solution.|'!|4|4|4|4First 
 5561  we observe that if 0|4|¬E|4|εn|4|¬W|4m,|'{A6}|h|εA|β(|βm|βα+
 5566  ↓|β1|β)|4|∂α=↓|4c|4α+↓|4|9A|βm|βn|4α+↓|4!!|42|gα_↓|gk|gα_↓|g
 5566  1A|β(|βm|βα_↓|βk|β)|βn|4α+↓|42|gα_↓|gmA|β0|βn|E|n|;
 5567  | A|β(|βm|βα+↓|β1|β)|βn|4|Lα=↓|4c|4α+↓|4|↔k|uc|)1|¬Ek|¬Em|)2
 5567  |gα_↓|gkA|β(|βm|βα+↓|β1|βα_↓|βk|β)|βn|4α+↓|42|gα_↓|gmA|β0|βn
 5567  >{A4}|Lα=↓|4c|4α+↓|4|f1|d32|)A|βm|βn|4α+↓|4|↔k|uc|)1|¬Ek|¬Wm
 5568  |)2|gα_↓|gk|gα_↓|g1A|β(|βm|βα_↓|βk|β)|βn|4α+↓|42|gα_↓|gmA|β0
 5568  |βn>{A4}|Lα=↓|4c|4α+↓|4|f1|d32|)A|βm|βn|4α+↓|4|f1|d32|)(A|βm
 5569  |βn|4α_↓|4c)>{A4}|Lα=↓|4|(c|d22|)|4α+↓|4A|βm|βn.>
 5571  {A6}|πHence |εA{U0}{H10L12M29}58320#Computer 
folio 432 galley 8
 5573  Programming!(A.-W./Knuth)!ch. 4!f. 432!g. 8|'
 5577  {A12}!|4|4|4|4Now let |εA|βm|4α=↓|4A|βm|βm. |πIf 
 5581  |εm|4|¬Q|40, |πwe have|'{A6}|εA|ur|)(mα+↓1)m|)|4|∂α=↓|4c|4α+
 5584  ↓|4|↔k|uc|)1|¬Ek|¬Emα+↓1|)|42|gα_↓|gkA|ur|)(mα+↓1α_↓k)m|)|4α
 5584  +↓|42|urα_↓mα_↓1|)|)A|β0|βm|'{A4}|Lα=↓|4c|4α+↓|4|f1|d32|)A|β
 5585  m|βm|4α+↓|4|↔k|uc|)1|¬Ek|¬Em|)|4{H12}({H10}2|urα_↓kα_↓1|)|)(
 5585  A|ur|)(mα_↓k)(mα+↓1)|)|4α_↓|4c/2){H12})|4{H10}α+↓|42|urα_↓mα
 5585  _↓1|)|)A|β0|βm>{A4}|Lα=↓|4c|4α+↓|4|f1|d32|)A|βm|βm|4α+↓|4|f1
 5586  |d32|)(A|ur|)(mα+↓1)(mα+↓1)|)|4α_↓|4a|4α_↓|42|gα_↓|gmb)|4α_↓
 5586  |4|f1|d34|)c(1|4α_↓|42|gα_↓|gm)>{A4}α+↓|42|urα_↓mα_↓1|)|)|↔a
 5587  |(c(m|4α+↓|41)|d22|)|4α+↓|4A|β0|β0|↔s|?{A4}|Lα=↓|4|f1|d32|)(
 5588  A|βm|4α+↓|4A|βm|βα+↓|β1)|4α+↓|4|f3|d34|)c|4α_↓|4|f1|d32|)a|4
 5588  α+↓|42|urα_↓mα_↓1|)|)(c|4α_↓|4b|4α+↓|4A|β0|β0)|4α+↓|4m2|gα_↓
 5588  |gm|gα_↓|g2c.>{A3}(36)|?{A6}|πSimilar maneuvering, 
 5592  as shown in exercise 19, establishes the relation|'
 5600  {A6}|εA|βn|βα+↓|β1|4α=↓|4|f3|d34|)A|βn|4α+↓|4|f1|d34|)A|βn|β
 5600  α_↓|β1|4α+↓|4|≤a|4α+↓|42|urα_↓nα_↓1|)|)|≤b|4α+↓|4(n|4α+↓|42)
 5600  2|urα_↓nα_↓1|)|)|≤g,!n|4|¬R|42,!(37)|?{A6}|πwhere|'
 5602  {A6}|≤a|4α=↓|4|f1|d34|)a|4α+↓|4|f7|d38|)c,|4|1|1|≤b|4α=↓|4A|
 5602  β0|β0|4α_↓|4b|4α_↓|4|f3|d32|)c,!!|πand!!|ε|≤g|4α=↓|4|f1|d32|
 5602  )c.|;{A6}|π!|4|4|4|4Thus the double recurrence 
 5607  (34) can be transformed into the single recurrence 
 5615  relation in (37). Use of the generating function 
 5623  |εG(z)|4α=↓|4A|β0|4α+↓|4A|β1z|4α+↓|4A|β2z|g2|4α+↓|1|1|¬O|4|¬
 5623  O|4|¬O |πnow transforms (37) into the equation|'
 5630  {A6}(1|4α_↓|4|f3|d34|)z|4α_↓|4|f1|d34|)z|g2)G(z)|4α=↓|4a|β0|
 5630  4α+↓|4a|β1z|4α+↓|4a|β2z|g2|4α+↓|4|(|≤a|d21|4α_↓|4z|)|4α+↓|4|
 5630  (|≤b|d21|4α_↓|4z/2|)|4α+↓|4|(|≤g|d2(1|4α_↓|4z/2)|g2|)|4,|;
 5631  {A3}(38)|?{A6}|πwhere |εa|β0, a|β1, |πand |εa|β2 
 5637  |πare constants that can be determined by the 
 5645  values of |εA|β0, A|β1, |πand |εA|β2. |π Since 
 5653  |ε1|4α_↓|4|f3|d34|)z|4α+↓|4|f1|d34|)z|g2|4α=↓|4(1|4α+↓|4|f1|
 5653  d34|)z)(1|4α_↓|4z), |πwe can express |εG(z) |πby 
 5659  the method of partial fractions in the form|'
 5667  {A6}|εG(z)|4α=↓|4b|β0|4α+↓|4b|β1z|4α+↓|4|(b|β2|d2(1|4α_↓|4z)
 5667  |g2|)|4α+↓|4|(b|β3|d21|4α_↓|4z|)|4α+↓|4|(b|β4|d2(1|4α_↓|4z/2
 5667  )|g2|)|4α+↓|4|(b|β5|d21|4α_↓|4z/2|)|4α+↓|4|(b|β6|d21|4α+↓|4z
 5667  /4|)|4.|;{A6}|πTedious manipulations produce 
 5671  the values of these constants |εb|β0,|4.|4.|4.|4,|4b|β6, 
 5677  |πand therefore the coe∃cients of |εG(z) |πare 
 5684  determined. We _nally obtain the solution|'{A6}|h|εA|βm|βn|4
 5690  |∂α=↓|4mc/2|4α+↓|4n(|9a|4α+↓|4|9c)|4α+↓|4(|f6|d325|)a|4α+↓|4
 5690  |9b|4α+↓|4|f7|d323|)c|4α+↓|4|9A|β0|β0)|4α+↓|42|gα_↓|gn(|9c)|
 5690  E|n|;| A|βn|βn|4|Lα=↓|4n(|f1|d35|)a|4α+↓|4|f7|d310|)c)|4α+↓|
 5691  4(|f16|d325|)a|4α+↓|4|f2|d35|)b|4α_↓|4|f23|d350|)c|4α+↓|4|f3
 5691  |d35|)A|β0|β0)>{A4}|L|4|4|4|4|4α+↓|42|gα_↓|gn(α_↓|f1|d33|)cn
 5692  |4α+↓|4|f2|d33|)b|4α_↓|4|f1|d39|)c|4α_↓|4|f2|d33|)A|β0|β0)>
 5693  {A4}|L|4|4|4|4|4α+↓|4(α_↓|f1|d34|))|gn(α_↓|f16|d325|)a|4α_↓|
 5693  4|f16|d315|)b|4α+↓|4|f16|d3225|)c|4α+↓|4|f16|d315|)A|β0|β0)|
 5693  4α+↓|4|f1|d32|)c|≤d|βn|β0;>{A6}| A|βm|βn|4|Lα=↓|4mc/2|4α+↓|4
 5694  n(|f1|d35|)a|4α+↓|4|f1|d35|)c)|4α+↓|4(|f6|d325|)a|4α+↓|4|f2|
 5694  d35|)b|4α+↓|4|f7|d350|)c|4α+↓|4|f3|d35|)A|β0|β0)|4α+↓|42|gα_
 5694  ↓|gn(|f1|d33|)c)>{A4}|L|4|4|4|4|4α+↓|4(α_↓|f1|d34|))|gn(α_↓|
 5695  f6|d325|)a|4α_↓|4|f2|d35|)b|4α+↓|4|f2|d375|)c|4α+↓|4|f2|d35|
 5695  )A|β0|β0),!!m|4|¬Q|4n.|J!(39)>{A6}|π!|4|4|4|4With 
 5697  these elaborate calculations done, we can readily 
 5704  determine the behavior of the lattice-point model. 
 5711  Assume that the inputs |εu |πand |εv |πto the 
 5720  algorithm are odd, and let |εm|4α=↓|4|"l|πlg|β2|4|εu|"L, 
 5726  n|4α=↓|4|"l|πlg|β2|4|εv|"L. |πThe average number 
 5730  of subtraction cycles, namely, the quantity |εC 
 5737  |πin the analysis of Program B, the number of 
 5746  times step B6 is executed, is obtained by setting 
 5755  |εa|4α=↓|41, b|4α=↓|40, c|4α=↓|41, A|β0|β0|4α=↓|41 
 5759  |πin the recurrence (34). By (39) we see that 
 5768  (for |εm|4|¬R|4n) |πthe _attice model predicts|'
 5774  {A6}C|4α=↓|4|f1|d32|)m|4α+↓|4|f2|d35|)n|4α+↓|4|f49|d350|)|4α
 5774  _↓|4|f1|d35|)|≤d|βm|βn|J!(40)|;{A6}|πsubtraction 
 5776  cycles, plus terms which rapidly go to zero as 
 5785  |εn |πapproaches in_nity.|'!|4|4|4|4The average 
 5790  number of times gcd(|εu,|4v)|4α=↓|41 |πis obtained 
 5796  by setting |εa|4α=↓|4b|4α=↓|4c|4α=↓|40, A|β0|β0|4α=↓|41; 
 5800  |πthis gives the probability that |εu |πand |εv 
 5808  |πare relatively prime, approximately |f3|d35|). 
 5813  Actually, since |εu |πand |εv |πare assumed to 
 5821  be odd, they should be relatively prime with 
 5829  probability |ε8/|≤p|g2 |π(see exercise 13), so 
 5835  this re⊗ects the degree of inaccuracy of the 
 5843  lattice-point model.|'!|4|4|4|4The average number 
 5848  of times a path from |ε(m,|4n) |πgoes through 
 5856  one of the ``diagonal'' points (|εm|¬S,|4m|¬S) 
 5862  |πfor some |εm|¬S|4|¬R|41 |πis obtained by setting 
 5869  |εa|4α=↓|41, b|4α=↓|4c|4α=↓|4A|β10*?*?*?4α=↓|40 
 5871  |πin (34); so we _nd this quantity is approximately|'
 5880  {A6}|ε|f1|d35|)n|4α+↓|4|f6|d325|)|4α+↓|4|f2|d35|)|≤d|βm|βn,!
 5880  !|πwhen!!|εm|4|¬R|4n.|;{A6}|πNow we can determine 
 5885  the average number of shifts, the number of times 
 5894  step B3 is performed. (This is the quantity |εD 
 5903  |πin Program B.) In any execution of Algorithm 
 5911  B, with |εu |πand |εv |πboth odd, the corresponding 
 5920  path in the lattice model must satisfy the relation|'
 5929  {A6}number|4of|4shifts|4α+↓|4number|4of|4diagonal|4points|4α
 5929  +↓|42|"llg|β2|4gcd(|εu,|4v)|"L|4α=↓|4m|4α+↓|4n,|;
 5930  {A6}|πsince we are assuming that |εr |πin (33) 
 5938  is always either 0 or 1. The average value of 
 5948  |"llg|β2|4gcd(|εu,|4v)|"L |πpredicted by the 
 5952  lattice-point model is approximately |f4|d35|) 
 5957  (see exercise 20). Therefore we have, for the 
 5965  total number of shifts,|'{A6}|εD|4|∂α=↓|4m|4α+↓|4n|4α_↓|4(|f
 5969  1|d35|)n|4α+↓|4|f6|d325|)|4α+↓|4|f2|d35|)|≤d|βm|βn)|4α_↓|f4|
 5969  d35|)|;{A4}|Lα=↓|4m|4α+↓|4|f4|d35|)n|4α_↓|4|f46|d325|)|4α_↓|
 5970  4|f2|d35|)|≤d|βm|βn,|J!(41)>{A6}|πwhen |εm|4|¬R|4n, 
 5973  |πplus terms which decrease rapidly to zero for 
 5981  large |εn.|'|π!|4|4|4|4To summarize the most 
 5987  important facts we have derived from the lattice-point 
 5995  model, we have shown that if |εu |πand |εv |πare 
 6005  odd and if |"llg|β2|4|εu|"L|4α=↓|4m, |"l|πlg|β2|4|εv|"L|4α=↓
 6009  |4n, |πthen the quantities |εC |πand |εD |πwhich 
 6017  are the critical factors in the running time 
 6025  of Program B will have average values given by|'
 6034  {A6}|εC|4α=↓|4|f1|d32|)m|4α+↓|4|f2|d35|)n|4α+↓|4O(1),!!D|4α=
 6034  ↓|4m|4α+↓|4|f4|d35|)n|4α+↓|4O(1),!!m|4|¬R|4n.|J!(42)|;
 6035  {A6}|πBut the model which we have used to derive 
 6044  (42) is only a pessimistic approximation to the 
 6052  true behavior; Table 1 compares the true average 
 6060  values of |εC, |πcomputed by actually running 
 6067  Algorithm B with all possible inputs, to the 
 6075  values predicted by the lattice-point model, 
 6081  for small |εm |πand |εn. |πThe lattice model 
 6089  is completely accurate when |εm |πor |εn |πis 
 6097  zero, but it tends to be less accurate when |≥≥|¬*?*?ends 
 6107  to be less accurate when |ε|¬Gm|4α_↓|4n|¬G |πis 
 6114  small and min(|εm,|4n) |πis large. When |εm|4α=↓|4n|4α=↓|49,
 6120   |πthe lattice-point model gives |εC|4α=↓|48.78, 
 6126  |πcompared to the true value |εC|4α=↓|47.58.|'
 6132  {A12}{A12}|π|∨T|∨a|∨b|∨l|∨e|9|∨1|;{A3}{H9L11M29}NUMBER|9OF|9
 6133  SUBTRACTIONS|9IN|9ALGORITHM|9B|;{A6}|∂!!|9|∂!!!!!!!!!!!!!!!|
 6134  ∂!!!|4|4|4|∂!!!!!!!!!!!!!!!|∂!!|9|∂|E|;|>|;|εn|;
 6138  |;n|;|;>{A6}|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!!|
 6142  4|4|4|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂!!|9|∂|E|;
 6143  |>|;0|;1|;2|;3|;4|;5|;|;0|;1|;2|;3|;4|;5|;|;>
 6160  |>0|;1.00|;2.00|;2.50|;3.00|;3.50|;4.00|;|;1.00|;
 6170  2.00|;2.50|;3.00|;3.50|;4.00|;0|;>|>1|;2.00|;
 6180  1.00|;2.50|;3.00|;3.50|;4.00|;|;2.00|;1.00|;3.00|;
 6189  2.75|;3.63|;3.94|;1|;>|>2|;2.50|;2.50|;2.25|;
 6199  3.38|;3.88|;4.38|;|;2.50|;3.00|;2.00|;3.50|;3.88|;
 6208  4.25|;2|;>|>3|;3.00|;3.00|;3.38|;3.25|;4.22|;
 6218  4.72|;|;3.00|;2.75|;3.50|;2.88|;4.13|;4.34|;3|;
 6227  >|>4|;3.50|;3.50|;3.88|;4.22|;4.25|;5.10|;|;3.50|;
 6238  3.63|;3.88|;4.13|;3.94|;4.80|;4|;>|>5|;4.00|;
 6248  4.00|;4.38|;4.72|;5.10|;5.19|;|;4.00|;3.94|;4.25|;
 6257  4.34|;4.80|;4.60|;5|;>{A3}|∂!!|9|∂!!!!!!!!!!!!!!!|∂!!!|4|4|4
 6262  |∂!!!!!!!!!!!!!!!|∂!!|9|∂|E|;|>m|;|πPredicted|4by|4model|;
 6266  |;Actual|4average|4values|;|εm|;>{A24}|H*?*?{U0}{H10L12M29}583
folio 434 galley 9
 6270  20#Computer Programming!(A.-W./Knuth)!ch. 4!f. 
 6273  434!g. 9|'{A12}|π!|4|4|4|4Empirical tests of 
 6278  Algorithm B with several million random inputs 
 6285  and with various values of |εm, n |πin the range 
 6295  29|4|¬E|4|εm,|4n|4|¬E|437 |πindicate that the 
 6299  actual average behavior of the algorithm is given 
 6307  by|'{A6}|h|εD|4|∂|¬V|4|9m|4α+↓|40.000n|4α+↓|40.0|4α⊗↓|40.7(0
 6308  .8)|gm|gα_↓|gn,!!m|4|¬R|4n.|E|n|; C|4|L|¬V|4|f1|d32|)m|4α+↓|
 6310  40.203n|4α+↓|41.9|4α_↓|40.4(0.6)|gm|gα_↓|gn,>
 6311  {A4}| D|4|L|¬V|4|9m|4α+↓|40.41n|9|1|4α_↓|40.5|4α_↓|40.7(0.6)
 6311  |gm|gα_↓|gn,!!m|4|¬R|4n.|J!(43)|;{A12}|πThese 
 6313  experiments showed a rather small standard deviation 
 6320  from the observed average values. The coe∃cients 
 6327  |f1|d32|) and 1 of |εm |πin (42) and (43) can 
 6337  be veri_ed rigorously without using the lattice-point 
 6344  approximation (see exercise 21), so the error 
 6351  in the lattice-point model is apparently in the 
 6359  coe∃cient of |εn |πwhich is too high.|'!|4|4|4|4The 
 6367  above calculations have been made under the assumption 
 6375  that |εu |πand |εv |πare odd and in the ranges 
 6385  |ε2|gm|4|¬W|4u|4|¬W|42|gm|gα+↓|g1, 2|gn|4|¬W|4v|4|¬W|42|gn|g
 6386  α+↓|g1. |πIf we say instead that |εu |πand |εv 
 6395  |πare to be |εany |πintegers, independently and 
 6402  uniformly distributed over the ranges|'{A6}|ε1|4|¬E|4u|4|¬W|
 6407  42|gN,!!1|4|¬E|4v|4|¬W|42|gN,|;{A6}|πthen we 
 6410  can calculate the average values of |εC |πand 
 6418  |εD |πfrom the data already given; in fact, if 
 6427  |εC|βm|βn |πdenotes the average value of |εC 
 6434  |πunder our earlier assumptions, exercise 22 
 6440  shows that we have|'{A6}|ε(2|gN|4α_↓|41)|g2C|4α=↓|4N|g2C|β0|
 6444  β0|4α+↓|42N|4|↔k|uc|)1|¬En|¬EN|)|4(N|4α_↓|4n)2|gn|gα_↓|g1C|β
 6444  n|β0|'{A4}α+↓|42|1|1|↔k|uc|)1|¬En|¬Wm|¬EN|)(N|4α_↓|4m)(N|4α_
 6445  ↓|4n)2|urmα+↓nα_↓2|)|)C|βm|βnα+↓|1|1|↔k|uc|)1|¬En|¬EN|)(N|4α
 6445  _↓|4n)|g22|ur2nα_↓2|)|)C|βn|βn,|J!(44)|'{A12}|πThe 
 6447  same formula holds for |εD |πin terms of |εD|βm|βn. 
 6456  |πIf the indicated sums are carried out using 
 6464  the approximations in (43), we obtain|'{A6}|εC|4|¬V|40.70N|4
 6470  α+↓|40(1),!!D|4|¬V|41.41N|4α+↓|4O(1).|;{A6}|π(See 
 6472  exercise 23.) This agrees perfectly with the 
 6479  results of further empirical tests, made on several 
 6487  million random inputs for |εN|4|¬E|430; |πthe 
 6493  latter tests show that we may take|'{A6}|εC|4α=↓|40.70N|4α_↓
 6500  |40.5,!!D|4α=↓|41.41N|4α_↓|42.7|J!(45)|;{A6}|πas 
 6502  good estimates of the values, given this distribution 
 6510  of the inputs |εu |πand |εv.|'|π!|4|4|4|4Richard 
 6517  Brent has discovered a continuous model which 
 6524  accounts for the leading terms in (45). Let us 
 6533  assume that |εu |πand |εv |πare large, and that 
 6542  the number of shifts per iteration has the value 
 6551  |εd |πwith probability exactly 2|ε|gd. |πIf we 
 6558  let |εX|4α=↓|4u/v, |πthe e=ect of steps B3<B5 
 6565  is to replace |εX |πby (|εX|4α_↓|41)/2|gd |πif 
 6572  |εX|4|¬Q|41, |πor by 2|ε|gd/(X|4α_↓|41) |πif 
 6577  |εX|4|¬W|41. |πThe random variable |εX |πhas 
 6583  a limiting distribution which makes it possible 
 6590  to analyze the average value of the ratio by 
 6599  which max(|εu,|4v) |πdecreases at each iteration; 
 6605  see exercise 25. Numerical calculations show 
 6611  that this maximum decreases by |ε|≤b|4α=↓|40.705971246 
 6617  |πbits per iteration; the agreement with experiment 
 6624  is so good that Brent's constant |ε|≤b |πmust 
 6632  be the true value of the number ``0.70'' in (45), 
 6642  and we should replace 0.203 by 0.206 in (43).|'
 6651  !|4|4|4|4This completes our analysis of the average 
 6658  values of |εC |πand |εD. |πThe other three quantities 
 6667  appearing in the running time of Algorithm B 
 6675  are rather easily analyzed; see exercises 6, 
 6682  7, and 8.|'!|4|4|4|4Thus we know approximately 
 6689  how Algorithm B behaves on the average. Let us 
 6698  now consider a ``worst case'' analysis: What 
 6705  values |"llg|β2|4|εu|"L|4α=↓|4m!!|πand!!|"llg|β2|4|εv|"L|4α=
 6706  ↓|4n,|;{A6}|πlet us try to _nd (|εu,|4v) |πwhich 
 6714  make the algorithm run most slowly. In view of 
 6723  the fact that the subtractions take somewhat 
 6730  longer than the shifts, when the auxiliary bookkeeping 
 6738  is considered, this question may be phrased by 
 6746  asking which |εu |πand |εv |πrequire most subtractions. 
 6754  The answer is somewhat surprising; the maximum 
 6761  value of |εC |πis exactly|'{A6}max(|εm,|4n)|4α+↓|41,|J!(46)|
 6766  ;{A6}|πalthough the lattice-point model would 
 6772  predict that substantially higher values of |εC 
 6779  |πare possible (see exercise 26). The derivation 
 6786  of the worst case (46) is quite interesting, 
 6794  so it has been left as an amusing problem for 
 6804  the reader to work out by himself (see exercises 
 6813  27, 28).|'{A24}|∨E|∨X|∨E|∨R|∨C|∨I|∨S|∨E|∨S|'{H9L11M29}{A12}|
 6816  9|1|≡1|≡.|9|4[|εM|*/|↔P|↔O|\] |πHow can (8), (9), 
 6821  (10), (11), and (12) be derived easily from (6) 
 6830  and (7)?|'{A3}|9|1|≡2|≡.|9|4[|εM|*/|↔P|↔P|\] |πGiven 
 6834  that |εu |πdivides |εv|β1v|β2|4.|4.|4.|4v|βn, 
 6838  |πprove that |εu |πdivides|'{A6}gcd(|εu,|4v|β1)|4|πgcd(|εu,|
 6842  4v|β2)|4|¬O|4|¬O|4|¬O|4|πgcd(|εu,|4v|βn).|;{A6}|π|9|1|≡3|≡.|
 6843  9|4|ε[M|*/|↔P|↔L|\] |πShow that the number of 
 6849  ordered pairs of positive integers |ε(u,|4v) 
 6855  |πsuch that lcm(|εu,|4v)|4α=↓|4n |πis the number 
 6861  of divisors of |εn|g2.|'{A3}|π|9|1|≡4|≡.|9|4|ε[M|*/|↔P|↔O|\] 
 6866  |πGiven positive integers |εu |πand |εv, |πshow 
 6873  that there are divisors |εu|¬S |πof |εu |πand 
 6881  |εv|¬S |πof |εv |πsuch that gcd(|εu|¬S,|4v|¬S)|4α=↓|41 
 6887  |πand |εu|¬Sv|¬S|4α=↓|4|πlcm(|εu,|4v).|'{A3}|π|9|1|≡5|≡.|9|4
 6889  |ε[M|*/|↔P|↔o|\] |πInvent an algorithm (analogous 
 6894  to algorithm B) for calculating the greatest 
 6901  common divisor of two integers based on their 
 6909  |εbalanced ternary |πrepresentation. Demonstrate 
 6913  your algorithm by applying it to the calculation 
 6921  of gcd(40902,24140).|'{A3}|9|1|≡6|≡.|9|4|ε[M|*/|↔P|↔P|\] 
 6924  |πGiven that |εu |πand |εv |πare random positive 
 6932  integers, _nd the mean and standard deviation 
 6939  of the quantity |εA (|πthe number of right shifts 
 6948  on both |εu |πand |εv |πduring the preparatory 
 6956  phase) which enters into the timing of Program 
 6964  B.|'{A3}|9|1|≡7|≡.|9|4|ε[M|*/|↔P|↔c|\] |πAnalyze 
 6967  the quantity |εB |πwhich enters into the timing 
 6975  of Program B.|'{A3}|9|1|≡8|≡.|9|4|ε[M|*/|↔P|↔M|\] 
 6979  |πShow that in Program B, the average value of 
 6988  |εE |πis approximately equal to |f1|d32|)|εC|π|βa|βv|βe, 
 6994  where |εC|π|βa|βv|βe is the average value of 
 7001  |εC.|'{A3}|π|9|1|≡9|≡.|9|4|ε[|*/|↔O|↔l|\] |πUsing 
 7004  Algorithm B and hand calculation, _nd gcd(31408,2718). 
 7011  Also _nd integers |εm |πand |εn |πsuch that 31408|εm|4α+↓|42
 7019  718n|4α=↓|4|πgcd(31408,2718), using Algorithm 
 7022  X.|'{A3}|≡1|≡0|≡.|9|4|ε[HM|*/|↔P|↔M|\] |πLet |εq|βn 
 7026  |πbe the number of ordered pairs of integers 
 7034  |ε(u,|4v) |πsuch that |ε1|4|¬E|4u,|4v|4|¬E|4n 
 7038  |πand gcd(|εu,|4v)|4α=↓|41. |πThe object of this 
 7044  exercise is to prove that lim|ε|βn|β|¬M|β|¬X|4q|βn/n|g2|4α=↓
 7049  |46/|≤p|g2, |πthereby establishing Theorem D.|'
 7054  !!|1|1a)|9Use the principle of inclusion and 
 7060  exclusion (Section 1.3.3) to show that|'{A6}|εq|βn|4α=↓|4n|g
 7066  2|4α_↓|4|↔k|uc|)p|β1|)|4|¬ln/p|β1|"L|g2|4α+↓|4|↔k|uc|)p|β1|¬
 7066  Wp|β2|)|4|"ln/p|β1p|β2|"L|g2|4α_↓|4|¬O|4|¬O|4|¬O|4,|;
 7067  {A6}|πwhere the sums are taken over all |εprime 
 7075  |πnumbers |εp|βi.|'|π!!|1|1b)|9The |εM|=4obius 
 7079  function |≤m(n) |πis de_ned by the rules |ε|≤m(1)|4α=↓|41, 
 7087  |≤m(p|β1p|β2|4.|4.|4.|4p|βr)|4α=↓|4(|→α_↓1)|gr 
 7088  |πif |εp|β1, p|β2,|4.|4.|4.|4,|4p|βr |πare distinct 
 7093  primes, and |ε|≤m(n)|4α=↓|40 |πif |εn |πis divisible 
 7100  by the square of a prime. Show that |εq|βn|4α=↓|4|¬K|βk|β|¬R
 7108  |β1|4|≤m(k)|"ln/k|"L|g2.|'|π!!|1|1c)|9As a consequence 
 7112  of (b), prove that lim|ε|βn|β|¬M|β|¬X|4q|βn/n|g2|4α=↓|4|¬K|β
 7116  k|β|¬R|β1|≤m(k)/k|g2.|'!!|1|1|πd)|9Prove that 
 7119  |ε(|¬K|βk|β|¬R|β1|≤m(k)/k|g2)(|¬K|βm|β|¬R|β1/m|g2)|4α=↓|41. 
 7120  Hint|*/: |\|πWhen the series are absolutely convergent 
 7127  we have|'{A6}|ε|↔a|↔k|uc|)k|¬R1|)|4a|βk/k|g2|↔s|↔a|↔k|uc|)m|
 7129  ¬R1|)|4b|βm/m|g2|↔s|4α=↓|4|↔k|uc|)n|¬R1|)|4|↔a|↔k|uc|)d|¬Dn|
 7129  )|4a|βdb|βn|β/|βd|↔s|↔hn|g2.|;{A6}|π|≡1|≡1|≡.|9|4|ε[M|*/|↔P|↔
 7130  P|\] |πWhat is the probability that gcd(|εu,|4v)|4|¬E|43? 
 7137  |π(See Theorem C.) What is the |εaverage |πvalue 
 7145  of gcd(|εu,|4v)?|'{A3}|π|≡1|≡2|≡.|9|4[|εM|*/|↔P|↔M|\] 
 7148  |π(E. Ces|=2aro.) If |εu |πand |εv |πare random 
 7156  positive integers, what is the average number 
 7163  of (positive) divisors they have in common? [|εHint|*/: 
 7171  |\|πSee the identity in exercise 10(d), with 
 7178  |εa|βk|4α=↓|4b|βm|4α=↓|41.]|'|H*?{U0}{H10L12M29}58320#Compute
folio 439 galley 10 WARNING: Some bad spots were on this tape.
 7179  r Programming!(A.-W./Knuth)!ch. 4!f. 439!g. 10|'
 7184  {H9L11M29}|π|≡1|≡3|≡.|9|4[|εHM|*/|↔P|↔L|\] |πGiven 
 7186  that |εu |πand |εv |πare random |εodd |πpositive 
 7194  integers, show that they are relatively prime 
 7201  with probability |ε8/|≤p|g2.|'{A3}|π|≡1|≡4|≡.|9|4|ε[M|*/|↔P|↔
 7204  o|\] |πWhat are the values of |εv|β1 |πand |εv|β2 
 7213  |πwhen algorithm X terminates?|'{A3}|≡1|≡5|≡.|9|4|ε[M|*/|↔P|↔
 7217  C|\] |πDesign an algorithm to |εdivide u by v 
 7226  modulo m, |πgiven positive integers |εu, v, |πand 
 7234  |εm, |πwith |εv |πrelatively prime to |εm. |πIn 
 7242  other words, your algorithm should _nd |εw, |πin 
 7250  the range |ε0|4|¬E|4w|4|¬W|4m, |πsuch that |εu|4|"o|4vw 
 7256  (|πmodulo |εm).|'{A3}|π|≡1|≡6|≡.|9|4[|*/|↔P|↔O|\] 
 7259  Use the text's method to _nd a general solution 
 7268  in integers to the following sets of equations:|'
 7276  {A6}!!|1|1a)|9|∂3|εx|4α+↓|47y|4α+↓|411z|4α=↓|41!!!!!!b)|9|∂3
 7276  x|4α+↓|47y|4α+↓|411z|4α=↓|4|4|4|4|41|'{A3}|L5x|4α+↓|47y|4α_↓
 7277  |4|9|15z|4α=↓|43|L5x|4α+↓|47y|4α_↓|4|9|15z|4α=↓|4|→α_↓3>
 7278  {A6}|π|≡1|≡7|≡.|9|4[|εM|*/|↔P|↔M|\] |πShow how 
 7281  Algorithm L can be extended (as Algorithm A was 
 7290  extended to Algorithm X) to obtain solutions 
 7297  of (15) when |εu |πand |εv |πare large.|'{A3}|≡1|≡8|≡.|9|4[|
 7305  εM|*/|↔L|↔p|\] |πLet |εu |πand |εv |πbe odd integers, 
 7313  independently and uniformly distributed in the 
 7319  ranges |ε2|gm|4|¬W|4u|4|¬W|42|gm|gα+↓|g1, 2|gn|4|¬W|4v|4|¬W|
 7321  42|gn|gα+↓|g1. |πWhat is the |εexact |πprobability 
 7327  that a single ``subtract and shift'' cycle in 
 7335  Algorithm B, namely, an operation that starts 
 7342  at step B6 and then stops after step B5 is _nished, 
 7353  reduces |εu |πand |εv |πto the ranges |ε2|gm|g|¬S|4|¬W|4u|4|
 7360  ¬W|42|gm|g|¬S|gα+↓|g1, 2|gn|g|¬S|4|¬W|4v|4|¬W|42|gn|g|¬S|gα+
 7361  ↓|g1, |πas a function of |εm, n, m|¬S, |πand 
 7370  |εn|¬S? |π(This exercise gives more accurate 
 7376  values for the transition probabilities than 
 7382  in the text's model.)|'{A3}|≡1|≡9|≡.|9|4[|εM|*/|↔P|↔M|\] 
 7387  |πComplete the text's derivation of (38) by establishing 
 7395  (37).|'{A3}|≡2|≡0|≡.|9|4|ε[M|*/|↔P|↔o|\] |πLet 
 7398  |ε|≤l|4α=↓|4|"l|πlg|4gcd(|εu,|4v)|"L. |πShow 
 7400  that the lattice-point model gives |ε|≤l|4α=↓|41 
 7406  |πwith probability |f1|d35|), |ε|≤l|4α=↓|42 |πwith 
 7411  probability |f1|d310|), |ε|≤l|4α=↓|43 |πwith 
 7415  probability |f1|d320|),|π etc., plus correction 
 7420  terms which go rapidly to zero as |εu, v |πapproach 
 7430  in_nity; hence the average value of |ε|≤l |πis 
 7438  approximately |f4|d35|). [|εHint|*/: |\|πconsider 
 7442  the relation between the probability of a path 
 7450  from |ε(m,|4n) |πto |ε(k|4α+↓|41, k|4α+↓|41) 
 7455  |πand a corresponding path from |ε(m|4α_↓|4k, 
 7461  n|4α_↓|4k) |πto (1,|41).]|'|≡2|≡1|≡.|9|4|ε[M|*/|↔P|↔C|\] 
 7465  |πLet |εC|βm|βn, D|βm|βn |πbe the average number 
 7472  of subtraction and shift cycles respectively, 
 7478  in Algorithm B when |εu |πand |εv |πare odd, 
 7487  |"llg|4|εu|"L|4α=↓|4m, |"l|πlg|4|εv|"L|4α=↓|4n. 
 7489  |πShow that for _xed |εn, C|βm|βn|4α=↓|4|f1|d32|)m|4α+↓|4O(1
 7494  ) |πand |εD|βm|βn|4α=↓|4m|4α+↓|4O(1) |πas |εm|4|¬M|4|¬X.|'
 7499  {A3}|π|≡2|≡2|≡.|9|4|ε[|*/|↔P|↔L|\] |πProve Eq. 
 7502  (44).|'{A3}|≡2|≡3|≡.|9|4|ε[M|*/|↔P|↔l|\] |πShow 
 7505  that if |εC|βm|βn|4α=↓|4|≤am|4α+↓|4|≤bn|4α+↓|4|≤g 
 7508  |πfor some constants |ε|≤a, |≤b, |πand |ε|≤g, 
 7515  |πthen|'{A6}|ε|↔k|uc|)1|¬En|¬Em|¬EN|)(N|4α_↓|4m)(N|4α_↓|4n)2
 7516  |gm|gα+↓|gn|gα_↓|g2C|βm|βn|4|∂α=↓|42|g2|gN{H12}({H9}|f11|d32
 7516  7|)(|≤a|4α+↓|4|≤b)N|4α+↓|4O(1){H12}){H10},|;{A6}| |↔k|uc|)1|
 7517  ¬En|¬EN|)(N|4α_↓|4n)|g22|g2|gn|gα_↓|g2C|βn|βn|4|Lα=↓|42|g2|g
 7517  N{H12}({H10}|f5|d327|)(|≤a|4α+↓|4|≤b)N|4α+↓|4O(1)}h12}){H10}
 7517  .>{A6}{H9L11M29}|π|≡2|≡4|≡.|9|4[|εM|*/|↔L|↔c|\] 
 7519  |πIf |εv|4α=↓|41 |πduring Algorithm B, while 
 7525  |εu |πis large, it may take fairly long for the 
 7535  algorithm to determine that gcd(|εu,|4v)|4α=↓|41. 
 7540  |πPerhaps it would be worth while to add a test 
 7550  at the beginning of step B5: ``If |εt|4α=↓|41, 
 7558  |πthe algorithm terminates with 2|ε|gk |πas the 
 7565  answer.'' Explore the question of whether {A3}|≡2|≡5|≡.|9|4|
 7571  ε[M|*/|↔P|↔C|\] |π(R. P. Brent.) Let |εu|βn, v|βn 
 7578  |πbe the values of |εu |πand |εv |πafter |εn 
 7587  |πiterations of steps B3<B5; let |εX|βn|4α=↓|4u|βn/v|βn, 
 7593  |πand assume that |εF|βn(x) |πis the probability 
 7600  that |εX|βn|4|¬E|4x, |πfor |ε0|4|¬E|4x|4|¬W|4|¬X. 
 7604  |π(a) Express |εF|βn|βα+↓|β1(x) |πin terms of 
 7610  |εF|βn(x), |πunder the assumption that step B4 
 7617  always branches to B3 with probability |f1|d32|). 
 7624  (b) Let |εG|βn(x)|4α=↓|4F|βn(x)|4α+↓|41|4α_↓|4F|βn(x|gα_↓|g1
 7626  ) |πbe the probability that |εY|βn|4|¬E|4x, |πfor 
 7633  |ε0|4|¬E|4x|4|¬E|41, |πwhere |εY|βn|4α=↓|4|πmin(|εu|βn,|4v|β
 7635  n)/|πmax(|εu|βn,|4v|βn). |πExpress |εG|βn|βα+↓|β1 
 7638  |πin terms of |εG|βn. |π(c) Express |εH|βn(x)|4α=↓|4|πPr(max
 7644  (|εu|βn|βα+↓|β1,|4v|βn|βα+↓|β1)/|πmax(|εu|βn,|4v|βn)|4|¬W|4x
 7644  ) |πin terms of |εG|βn.|'{A3}|π|≡2|≡6|≡.|9|4|ε[M|*/|↔P|↔L|\] 
 7650  |πWhat is the length of the longest path from 
 7659  |ε(m,|4n) |πto (0,|40) in the lattice-point model?|'
 7666  {A3}|≡2|≡7|≡.|9|4|ε[M|*/|↔P|↔l|\] |πFind values 
 7669  of |εu, v |πwith |"llg|4|εu|"L|4α=↓|4|εm, |"l|πlg|4|εv|"L|4α
 7674  =↓|4|εn, |πsuch that Algorithm B requires |εm|4α+↓|41 
 7681  |πsubtraction steps given |εm|4|¬R|4n|4|¬R|41.|'
 7685  {A|≡3|≡3|≡.|9|4|ε[M|*/|↔M|↔o|\] |πAnalyze V. C. 
 7689  Harris's ``binary Euclidean algorithm.''|'{A3}|≡3|≡4|≡.|9|4|
 7693  ε[M|*/|↔L|↔P|\] |π(R. W. Gosper.) Demonstrate 
 7698  how to modify Algorithm B for large numbers using 
 7707  ideas analogous to those in Algorithm L.|'{A18}{H10L12M29}{J
 7714  1}α/↓|∨4|∨.|∨5|∨.|∨3|∨.|9|4|∨A|∨n|∨a|∨l|∨y|∨s|∨i|∨s|9|∨o|∨f|
 7714  9|∨E|∨u|∨c|∨l|∨i|∨d|∨'|∨s|9|∨A|∨l|∨g|∨o|∨r|∨i|∨t|∨h|∨m|'
 7715  {A6}The execution time of Euclid's algorithm 
 7721  depends on |εT, |πthe number of times the division 
 7730  step A2 is performed. (See Algorithm 4.5.2A and 
 7738  Program 4.5.2A.) The same quantity |εT |πis an 
 7746  important factor in the running time of other 
 7754  algorithms, for example, the evaluation of functions 
 7761  satisfying a reciprocity formula (see Section 
 7767  3.3.3). We shall see in this section that the 
 7776  mathematical analysis of this quantity |εT |πis 
 7783  interesting and instructive.|'{A12}|≡R|≡e|≡l|≡a|≡t|≡i|≡o|≡n 
 7787  |≡t|≡o |≡c|≡o|≡n|≡t|≡i|≡n|≡u|≡e|≡d |≡f|≡r|≡a|≡c|≡t|≡i|≡o|≡n|
 7789  ≡s|≡.|9|4Euclid's algorithm is intimately connected 
 7794  with |εcontinued fractions, |πwhich are expressions 
 7800  of the form|'{A6}|h|εa|β1|4α+↓|4|∂a|β2|4α+↓|4b|β3|∂.|4.|4.|∂
 7803  a|βn|βα_↓|β1|4α+↓|4b|βn|4|∂α=↓{A12}a|β1|4α+↓|4b|β2|'
 7804  |La|β2|4α+↓|4b|β3|'|L|L|L|Lα=↓|4b|β1/{H12}({H10}a|β1|4α+↓|4b
 7805  |β2/(a|β2|4α+↓|4b|β3/(|¬O|4|¬O|4|¬O|4/(a|βn|βα_↓|β1|4α+↓|4b|
 7805  βn/a|βn)|4|¬O|4|¬O|4|¬O)){H12}){H10}.>|L|L|¬O|4|¬O|4|¬O>
 7807  {A12}|L|L|La|βn|βα_↓|β1|4α+↓|4b|βn>{A7}|L|L|L|La|βn|J!(1)>
 7809  |Hom{U0}{H10L12M29}58320#Computer Programming!(A.-W./Knuth)!
folio 441 galley 11
 7810  ch. 4!f. 441!g. 11|'{A12}|πContinued fractions 
 7816  have a beautiful theory which is the subject 
 7824  of several books, e.g., O. Perron, |εDie Lehre 
 7832  von den Kettenbr|=4uchen, |π3rd ed. (Stuttgart: 
 7838  Teubner, 1954), 2 vols.; A. Khinchin, |εContinued 
 7845  fractions, |πtr. by Peter Wynn (Groningen: P. 
 7852  Noordho=, 1963); H. S. Wall, |εAnalytic theory 
 7859  of continued fractions (|πNew York: Van Nostrand, 
 7866  1948); see also J. Tropfke, |εGeschichte der 
 7873  Elementar-Mathematik |≡6 (|πBerlin: Gruyter, 
 7877  1924), 74<84, for the early history of the subject. 
 7886  It is necessary to limit ourselves to a comparatively 
 7895  brief treatment of the theory here, studying 
 7902  only those aspects which give us more insight 
 7910  into the behavior of Euclid's algorithm.|'!|4|4|4|4The 
 7917  continued fractions which are of primary interest 
 7924  to us are those in which all the |εb'|πs in (1) 
 7935  are equal to unity. For convenience in notation, 
 7943  let us de_ne|'{A6}|"C|εx|β1,|4x|β2,|4.|4.|4.|4,|4x|βn|"C|4α=
 7946  ↓|41/{H12}({H10}x|β1|4α+↓|41/(x|β2|4α+↓|41/(|¬O|4|¬O|4|¬O|4α
 7946  +↓|41/(x|βn)|4|¬O|4|¬O|4|¬O)){H12}){H10}.|J!(2)|;
 7947  {A6}|πIf |εn|4α=↓|40, |πthe symbol |"C|εx|β1,|4.|4.|4.|4,|4x
 7951  |βn|"C |πis taken to mean 0. Thus, for example,|'
 7960  {A6}|"C|εx|β1|"C|4α=↓|4|(1|d2x|β1|)|4,!!|"Cx|β1,|4x|β2|"C|4α
 7960  =↓|4|(1|d2x|β1|4α+↓|4(1/x|β2)|)|4α=↓|4|(x|β2|d2x|β1x|β2|4α+↓
 7960  |41|)|4.|J!(3)|;{A6}|πLet us also de_ne the polynomials 
 7967  |εQ|βn(x|β1,|4x|β2,|4.|4.|4.|4,|4x|βn) |πof |εn 
 7970  |πvariables, for |εn|4|¬R|40, |πby the rule|'
 7976  {A6}|ε|h|εQ|βn(x|β1,|4x|β2,|4.|4.|4.|4,|4x|βn)|4α=↓|4!|∂juio
 7976  khyujnmk,|E|n|'|L1,|J!|πif!|4|4|εn|4α=↓|40;>{A4}Q|βn(x|β1,|4
 7978  x|β2,|4.|4.|4.|4,|4x|βn)|4α=↓|4!x|β1,|J!|πif!|4|4|εn|4α=↓|41
 7978  ;|'{A4}|Lx|β1Q|βn|βα_↓|β1(x|β2,|4.|4.|4.|4,|4x|βn)|4α+↓|4Q|β
 7979  n|βα_↓|β2(x|β3,|4.|4.|4.|4,|4x|βn),|J!|πif!|9|εn|4|¬Q|41.>
 7980  {A3}(4)|?|π{A6}Thus |εQ|β2(x|β1,|4x|β2)|4α=↓|4x|β1x|β2|4α+↓|
 7982  41, Q|β3(x|β1,|4x|β2,|4x|β3)|4α=↓|4x|β1x|β2x|β3|4α+↓|4x|β1|4
 7983  α+↓|4x|β3, |πetc. In general, as noted by L. 
 7991  Euler in the wighteenth century, |εQ|βn(x|β1,|4x|β2,|4.|4.|4
 7996  .|4,|4x|βn) |πis the sum of all terms obtainable 
 8004  by starting with |εx|β1x|β2|4.|4.|4.|4x|βn |πand 
 8009  deleting zero or more nonoverlapping pairs of 
 8016  consecutive variables |εx|βjx|βj|βα+↓|β1; |πthere 
 8020  are |εF|βn|βα+↓|β1 |πsuch terms. The polynomials 
 8026  de_ned in (4) are called ``continuants.``|'!|4|4|4|4The 
 8033  basic property of the |εQ-|πpolynomials is that|'
 8040  {A6}|"C|εx|β1,|4x|β2,|4.|4.|4.|4,|4x|βn|"C|4α=↓|4Q|βn|βα_↓|β
 8040  1(x|β2,|4.|4.|4.|4,|4x|βn)/Q|βn(x|β1,|4x|β2,|4.|4.|4.|4,|4x|
 8040  βn),!!n|4|¬R|41.|J!(5)|;{A6}|πThis can be proved 
 8045  by induction, since it implies that|'{A6}|εx|β0|4α+↓|4|"Cx|β
 8051  1,|4.|4.|4.|4,|4x|βn|"C|4α=↓|4Q|βn|βα+↓|β1(x|β0,|4x|β1,|4.|4
 8051  .|4.|4,|4x|βn)/Q|βn(x|β1,|4.|4.|4.|4,|4x|βn);|;
 8052  {A6}|πhence |ε|"Cx|β0,|4x|β1,|4.|4.|4.|4,|4x|βn|"C 
 8054  |πis the reciprocal of the latter quantity.|'
 8061  !|4|4|4The |εQ-|πpolynomials are symmetrical 
 8065  in the sense that|'{A6}|εQ|βn(x|β1,|4x|β2,|4.|4.|4.|4,|4x|βn
 8069  )|4α=↓|4Q|βn(x|βn,|4.|4.|4.|4,|4x|β2,|4x|β1).|J!(6)|;
 8070  {A6}|πThis follows from Euler's observation above, 
 8076  and as a consequence we have|'{A6}|εQ|βn(x|β1,|4.|4.|4.|4,|4
 8082  x|βn)|4α=↓|4x|βnQ|βn|βα_↓|β1(x|β1,|4.|4.|4.|4,|4x|βn|βα_↓|β1
 8082  )|4α+↓|4Q|βn|βα_↓|β2(x|β1,|4.|4.|4.|4,|4x|βn|βα_↓|β2).|J!(7)
 8082  |;|π{A6}The |εQ-|πpolynomials also satisfy the 
 8088  important identity|'{A6}|ε!Q|βN(x|β1,|4.|4.|4.|4,|4x|βn)Q|βn
 8090  (x|β2,|4.|4.|4.|4,|4x|βn|βα+↓|β1)|4α_↓|4Q|βn|βα+↓|β1(x|β1,|4
 8090  .|4.|4.|4,|4x|βn|βα+↓|β1)Q|βn|βα_↓|β1(x|β2,|4.|4.|4.|4,|4x|β
 8090  n)|'{A4}α=↓|4(|→α_↓1)|gn,!!n|4|¬R|41.!(8)|?{A6}|π(See 
 8093  exercise 4.) The latter equation in connection 
 8100  with (5) implies that|'{A6}!|ε|"Cx|β1,|4.|4.|4.|4,|4x|βn|"C|
 8104  4α=↓|4|(1|d2q|β0q|β1|)|4α_↓|4|(1|d2q|β1q|β2|)|4α+↓|4|(1|d2q|
 8104  β2q|β3|)|4α_↓|4|¬O|4|¬O|4|¬O|4α+↓|4|((|→α_↓1)|gn|gα_↓|g1|d2q
 8104  |βn|βα_↓|β1q|βn|)|4,|'|πwhere!!|εq|βk|4α=↓|4Q|βk(x|β1,|4.|4.
 8105  |4.|4,|4x|βk).!(9)|?{A6}|πThus the |εQ-|πpolynomials 
 8109  are intimately related to continued fractions.|'
 8115  !|4|4|4|4Every real number |εX |πin the range 
 8122  |ε0|4|¬E|4X|4|¬W|41 |πhas a |εregular continued 
 8127  fraction |πde_ned as follows: Let |εX|β0|4α=↓|4X, 
 8133  |πand for all |εn|4|¬R|40 |πsuch that |εX|βn|4|=|↔6α=↓|40 
 8140  |πlet|'{A6}|εA|βn|βα+↓|β1|4α=↓|4|"l1/X|βn|"L,!!X|βn|βα+↓|β1|
 8141  4α=↓|41/X|βn|4α_↓|4A|βn|βα+↓|β1.|J!(10)|;{A6}|πIf 
 8143  |εX|βn|4α=↓|40, |πthe quantities |εA|βn|βα+↓|β1 
 8147  |πand |εX|βn|βα+↓|β1 |πare not de_ned, and the 
 8154  regular continued fraction for |εX |πis |ε|"CA|β1,|4.|4.|4.|
 8160  4,|4A|βn|"C. |πIf |εX|4|=|↔6α=↓|40, |πthis de_nition 
 8165  guarantees that |ε0|4|¬E|4X|βn|βα+↓|β1|4|¬W|41, 
 8168  |πso each of the |εA'|πs is a positive integer. 
 8177  The de_nition (10) clearly implies that|'|ε{A6}X|4α=↓|4X|β0|
 8183  4α=↓|4|(1|d2A|β1|4α+↓|4X|β1|)|4α=↓|4|(1|d2A|β1|4α+↓|41/(A|β2
 8183  |4α+↓|4X|β2)|)|4α=↓|4|¬O|4|¬O|4|¬O|4;|;{A6}|πhence|'
 8185  {A6}|εX|4α=↓|4|"CA|β1,|4.|4.|4.|4,|4A|βn|βα_↓|β1,|4A|βn|4α+↓
 8185  |4X|βn|"C|J!(11)|;{A6}|πfor all |εn|4|¬R|41, 
 8189  |πwhenever |εX|βn |πis de_ned. Therefore, if 
 8195  |εX|βn|4α=↓|40, X|4α=↓|4|"CA|β1,|4.|4.|4.|4,|4A|βn|"C. 
 8197  |πIf |εX|βn|4|=|↔6α=↓|40, X |πlies |εbetween 
 8202  |"CA|β1,|4.|4.|4.|4,|4A|βn|"C |πand |ε|"CA|β1,|4.|4.|4.|4,|4
 8204  A|βn|4α+↓|41|"C, |πsince by (7) the quantity 
 8210  |εq|βn|4α=↓|4Q|βn(A|β1,|4.|4.|4.|4,|4A|βn|4α+↓|4X|βn) 
 8211  |πincreases monotonically from |εQ|βn(A|β1,|4.|4.|4.|4,|4A|β
 8214  n) |πup to |εQ|βN(A|β1,|4.|4.|4.|4,|4A|βn|4α+↓|41) 
 8218  |πas |εX|βn |πincreases from 0 to 1, and by (9) 
 8228  the continued fraction increases or decreases 
 8234  when |εq|βn |πincreases, according as |εn |πis 
 8241  even or odd.. In fact,|'{A6}{H10L12M29}|¬G|εX|4α_↓|4|"CA|β1,
 8246  |4.|4.|4.|4,|4A|βn|"C|¬G|4|∂α=↓|4|¬G|"CA|β1,|4.|4.|4.|4,|4A|
 8246  βn|4α+↓|4X|βn|"C|4α_↓|4|"CA|β1,|4.|4.|4.|4,|4A|βn|"C|¬G|'
 8247  {A3}|Lα=↓|4|¬G|"CA|β1,|4.|4.|4.|4,|4A|βn,|41/X|βn|"C|4α_↓|4|
 8247  "CA|β1,|4.|4.|4.|4,|4A|βn|"C|¬G>{A3}|Lα=↓|4|↔g|(Q|βn(A|β2,|4
 8248  .|4.|4.|4,|4A|βn,|41/X|βn)|d2Q|βn|βα+↓|β1(A|β1,|4.|4.|4.|4,|
 8248  4A|βn,|41/X|βn)|)|4α_↓|4|(Q|βn|βα_↓|β1(A|β2,|4.|4.|4.|4,|4A|
 8248  βn|d2Q|βn(A|β1,|4.|4.|4.|4,|4A|βn)|)|1|↔g>{A3}|Lα=↓|41/Q|βn(
 8249  A|β1,|4.|4.|4.|4,|4A|βn)Q|βn|βα+↓|β1(A|β1,|4.|4.|4.|4,|4A|βn
 8249  ,|41/X|βn)>{A3}|L|¬E|41/Q|βn(A|β1,|4.|4.|4.|4,|4A|βn)Q|βn|βα
 8250  +↓|β1(A|β1,|4.|4.|4.|4,|4A|βn,|4A|βn|βα+↓|β1)|J!(12)>
 8251  {A6}|πby (5), (8), and (10). Therefore |"C|εA|β1,|4.|4.|4.|4
 8257  ,|4A|βn|"C |πis an extremely close approximation 
 8263  to |εX. |πIf |εX |πis irrational, it is impossible 
 8272  to have |εX|βn|4α=↓|40 |πfor any |εn, |πso the 
 8280  regular continued fraction expansion in this 
 8286  case is an |εin⊂nite continued fraction |"CA|β1,|4A|β2,|4A|β
 8292  3,|4.|4.|4.|"C. |πThe value of such a continued 
 8299  fraction is de_ned to be lim|ε|βn|β|¬M|β|¬X|"CA|β1,|4A|β2,|4
 8304  .|4.|4.|4,|4A|βn|"C, |πand from the inequality 
 8309  (12) it is clear that this limit equals |εX.|'
 8318  |π!|4|4|4|4The regular continued fraction expansion 
 8323  of real numbers has several properties analogous 
 8330  to the representation of numbers in the decimal 
 8338  system. If we use the formulas above to compute 
 8347  the regular continued fraction expansions of 
 8353  some familiar real numbers, we _nd, for example, 
 8361  that|'{A6}{H9L11M29}|h|¬H|f8|d329|4|∂α=↓|4kojiuh|E|n|'
 8363  | |f8|d329|)|4|Lα=↓|4|"C3,|31,|31,|31,|32|"C;>
 8364  {B3}| {U19}|f8|d329|)|)|4|Lα=↓|4|"C1,|31,|39,|32,|32,|33,|32
 8364  ,|32,|39,|31,|32,|31,|39,|32,|32,|33,|32,|32,|39,|31,|32,|31
 8364  ,|39,|32,|32,|33,|32,|32,|39,|31,|3.|3.3.|"C;|'
 8365  {A2}| |=|g3|¬H|v42|)|4|Lα=↓|41|4α+↓|4|"C3,|31,|35,|31,|31,|3
 8365  4,|31,|31,|38,|31,|314,|31,|310,|32,|31,|34,|312,|32,|33,|32
 8365  ,|31,|33,|34,|31,|31,|32,|314,|33,|3.|3.|3.|"C;>
 8366  {A2}|ε| |≤p|4|Lα=↓|43|4α+↓|4|"C7,|315,|31,|3292,|31,|31,|31,
 8366  |32,|31,|33,|31,|31,|314,|32,|31,|31,|32,|32,|32,|32,|31,|38
 8366  4,|32,|32,|31,|31,|315,|33,|313,|3.|3.|3.|"C;>
 8367  {A2}| e|4|Lα=↓|42|4α+↓|4|"C1,|32,|31,|31,|34,|31,|31,|36,|31
 8367  ,|31,|38,|31,|31,|312,|31,|31,|312,|31,|31,|314,|31,|31,|316
 8367  ,|31,|31,|318,|31,|3.|3.|3.|"C;>| |≤g|4|Lα=↓|4|"C1,|31,|32,|
 8368  31,|32,|31,|34,|33,|313,|35,|31,|31,|38,|31,|32,|34,|31,|31,
 8368  |340,|31,|311,|33,|37,|31,|37,|31,|31,|35,|3.|3.|3.|"C;>
 8369  {A2}| |≤f|4|Lα=↓|41|4α+↓|4|"C1,|31,|31,|31,|31,|31,|31,|31,|
 8369  31,|3.|3.|3.|"C.|J!(13)>{A6}{H10L12M29}|πThe 
 8371  numbers |εA|β1,|4A|β2,|4.|4.|4. |πare called 
 8375  the |εpartial quotients |πof |εX. |πNote the 
 8382  regular pattern which appears in the partial 
 8389  quotients for |¬H|v48/29|), |ε|≤f, |πand |εe; 
 8395  |πthe reasons for this behavior are discussed 
 8402  in exercises 12 and 16. There is no apparent 
 8411  pattern in the partial quotients for |=|g3|¬H|v42|), 
 8418  |ε|≤p, |πor |ε|≤g.|'|π!|4|4|4|4It is interesting 
 8424  to note that the ancient Greek's _rst de_nition 
 8432  of real numbers, once they had discovered the 
 8440  existence of irrationals, was essentially in 
 8446  terms of in_nite continued fractions. (Later 
 8452  they adopted the suggestion of Eudoxus that |εx|4α=↓|4y 
 8460  |πshould be de_ned instead as ``|εx|4|¬W|4r |πif 
 8467  and only if |εy|4|¬W|4r, |πfor all rational |εr.'' 
 8475  |πSee O. Becker, |εQuellen und Studien zur Geschichte 
 8483  Math., Astron., Physik |π(B) |≡2 (1933), 311<333.)|'
folio 444 galley 12
 8490  |Hβ*?*?*?*?{U0}{H10L12M29}58320#Computer Programming!(A.-W./Knut
 8491  h)!ch. 4!f. 444!g. 12|'{A12}!|4|4|4|4When |εX 
 8497  |πis a rational number, the regular continued 
 8504  fraction corresponds in a natural way to Euclid's 
 8512  algorithm. Let us assume that |εX|4α=↓|4v/u, 
 8518  |πwhere |εu|4|¬Q|4v|4|¬R|40. |πThe regular continued 
 8523  fraction process starts with |εX|β0|4α=↓|4X; 
 8528  |πlet us de_ne |εU|β0|4α=↓|4u, V|β0|4α=↓|4v. 
 8533  |πAssuming that |εX|βn|4α=↓|4V|βn/U|βn|4|=|↔6α=↓|40, 
 8536  (10) |πbecomes |'{A6}|εA|βn|βα+↓|β1|4α=↓|4|"lU|βn/V|βn|"L,|;
 8540  {B7}(14)|?{A11}X|βn|βα+↓|β1|4α=↓|4U|βn/V|βn|4α_↓|4A|βn|βα+↓|
 8541  β1|4α=↓|4(U|βn|4|πmod|4|εV|βn)/V|βn.|;{A6}|πTherefore, 
 8543  if we de_ne|'{A6}|εU|βn|βα+↓|β1|4α=↓|4V|βn,!!V|βn|βα+↓|β1|4α
 8546  =↓|4U|βn|4|πmod|4|εV|βn,|J!(15)|;{A6}|πthe condition 
 8549  |εX|βn|4α=↓|4V|βn/U|βn |πholds throughout the 
 8553  process. Furthermore, (15) is precisely the transformation 
 8560  made on the variables |εu, v |πin Euclid's algorithm 
 8569  (see Algorithm 4.5.2A, step A2). For example, 
 8576  since |f8|d329|)|4α=↓|4|"C3, 1, 1, 2|"C, we know 
 8583  that Euclid's algorithm applied to |εu|4α=↓|429, 
 8589  v|4α=↓|48 |πwill require exactly _ve division 
 8595  steps, and the quotients |ε|"lu/v|"L |πin step 
 8602  A2 will be successively 3, 1, 1, 1, and 2. Note 
 8613  that the last partial quotient |εA|βn |πmust 
 8620  be 2 or more when |εX|βn|4α=↓|40, n|4|¬R|41, 
 8627  |πsince |εX|βn|βα_↓|β1 |πis less than unity.|'
 8633  !|4|4|4|4From this correspondence with Euclid's 
 8638  algorithm we can see that the regular continued 
 8646  fraction for |εX |πterminates at some step with 
 8654  |εX|βn|4α=↓|40 |πif and only if |εX |πis rational; 
 8662  for it is obvious that |εX|βn |πcannot be zero 
 8671  if |εX |πis irrational, and, conversely, we know 
 8679  that Euclid's algorithm always terminates. If 
 8685  the partial quotients obtained during Euclid's 
 8691  algorithm are |εA|β1,|4A|β2,|4.|4.|4.|4,|4A|βn, 
 8694  |πthen we have, by (5),|'{A6}|ε|(v|d2u|)|4α=↓|4|(Q|βn|βα_↓|β
 8699  1(A|β2,|4.|4.|4.|4,|4A|βn)|d2Q|βn(A|β1,|4A|β2,|4.|4.|4.|4,|4
 8699  A|βn)|)|4.|J!(16)|;{A6}|πThis formula holds also 
 8704  if Euclid's algorithm is applied for |εu|4|¬E|4v, 
 8711  |πwhen |εA|β1|4α=↓|40. |πFurthermore, because 
 8715  of (8), |εQ|βn|βα_↓|β1(A|β2,|4.|4.|4.|4,|4A|βn) 
 8718  |πand |εQ|βn(A|β1,|4A|β2,|4.|4.|4.|4,|4A|βn) 
 8720  |πare relatively prime, and the fraction on the 
 8728  right-hand side of (16) is in lowest terms; therefore|'
 8737  {A6}|εu|4α=↓|4Q|βn(A|β1,|4A|β2,|4.|4.|4.|4,|4A|βn)d,!!v|4α=↓
 8737  |4Q|βn|βα_↓|β1(A|β2,|4.|4.|4.|4,|4A|βn)d,|J!(17)|;
 8738  {A6}|πwhere |εd|4α=↓|4|πgcd(|εu,|4v).|'{A12}|π|≡T|≡h|≡e 
 8741  |≡w|≡o|≡r|≡s|≡t |≡c|≡a|≡s|≡e|≡.|9|4|πWe can now 
 8745  apply these observations to determine the behavior 
 8752  of Euclid's algorithm in the ``worst case,'' 
 8759  or in other|words to give an upper bound on the 
 8769  number of division steps. The worst case occurs 
 8777  when the inputs are consecutive Fibonacci numbers:|'
 8784  {A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡F (G. Lam|=1e, 1845). 
 8789  |εFor n|4|¬R|41, let u and v be integers with 
 8798  u|4|¬Q|4v|4|¬Q|40 such that Euclid's algorithm 
 8803  applied to u and v requires exactly n division 
 8812  steps, and such that u js as small as possible 
 8822  satisfying these conditions. Then u|4α=↓|4F|βn|βα+↓|β2, 
 8827  v|4α=↓|4F|βn|βα+↓|β1.|'{A6}Proof.|9|4|πBy (17), 
 8830  we must have |εu|4α=↓|4Q|βn(A|β1,|4A|β2,|4.|4.|4.|4,|4A|βn)d
 8833   |πwhere |εA|β1, A|β2, . . . , A|βn, d |πare 
 8844  positive integers and |εA|βn|4|¬R|42. |πSince 
 8849  |εQ|βn |πis a polynomial with nonnegative coe∃cients, 
 8856  involving all of the variables, the minimum value 
 8864  is values in (17) yields the desired result.|'
 8872  {A12}!|4|4|4|4(This theorem has the historical 
 8877  claim of being the _rst prictical application 
 8884  of the Fibonacci sequence; since then many other 
 8892  applications of the Fibonacci numbers to algorithms 
 8899  and to the study of algorithms have been discovered.)|'
 8908  !|4|4|4|4As a consequence of Theorem F we have 
 8916  an important corollary:|'{A12}|≡C|≡o|≡r|≡o|≡l|≡l|≡a|≡r|≡y|≡.
 8919  |9|4|εIf 0|4|¬E|4u, v|4|¬W|4N, the number of 
 8925  division steps required when Algorithm 4.5.2|πA 
 8931  |εis applied to u, v is at most |"plog|β|≤f|4(|¬H|v45|)N)|"P
 8939  |4α_↓|42.|'{A12}Proof.|9|4|πBy Theorem F, the 
 8944  maximum number of steps, |εn, |πoccurs when |εu|4α=↓|4F|βn 
 8952  |πand |εv|4α=↓|4F|βn|βα+↓|β1, |πwhere |εn |πis 
 8957  as large as possible with |εF|βn|βα+↓|β1|4|¬W|4N. 
 8963  (|πThe _rst division step in this case merely 
 8971  interchanges |εu |πand |εv |πwhen |εn|4|¬Q|41.) 
 8977  |πSince |εF|βn|βα+↓|β1|4|¬W|4N, |πwe have |ε|≤f|gn|gα+↓|g1/|
 8981  ¬H|v45|)|4|¬W|4N (|πsee Eq. 1.2.8<15), so |εn|4α+↓|41|4|¬W|4
 8986  |πlog|ε|β|≤f|4(|¬H|v45|)N).|'{A12}|πNote that 
 8989  log|ε|β|≤f|4(|¬H|v45|)N) |πis approximately 2.078|4ln|4|εN|4
 8992  α+↓|41.672|4|¬V|44.785|4|πlog|β1|β0|4|εN|4α+↓|41.672. 
 8993  |πSee exercises 31 and 36 for extensions of Theorem 
 9002  F.|'{A12}|≡A|≡n |≡a|≡p|≡p|≡r|≡o|≡x|≡i|≡m|≡a|≡t|≡e 
 9005  |≡m|≡o|≡d|≡e|≡l|≡.|9|4Now that we know the maximum 
 9011  number of division steps that can occur, let 
 9019  us attempt to _nd the |εaverage |πnumber. Let 
 9027  |εT(m,|4n) |πbe the number of division steps 
 9034  that occur when |εu|4α=↓|4m, v|4α=↓|4n |πare 
 9040  input to Euclid's algorithm. Thus|'{A6}|εT(m,|40)|4α=↓|40;!!
 9045  T(m,|4n)|4α=↓|41|4α+↓|4T(n,|4m|4|πmod|4|εn)!!|πif!!|εn|4|¬R|
 9045  41.|J!(18)|;{A6}|πLet |εT|βn |πbe the average 
 9051  number of division steps when |εv|4α=↓|4n |πand 
 9058  when |εu |πis chosen at random; since only the 
 9067  value of |εu |πmod |εv |πa=ects the algorithm 
 9075  after the _rst division step, we may write|'{A6}|εT|βn|4α=↓|
 9083  4|(1|d2n|)|4|↔k|uc|)0|¬Ek|¬Wn|)|4T(k,|4n).|J!(19)|;
 9084  {A6}|πFor example, |εT(0,|45)|4α=↓|41, T(1,|45)|4α=↓|42, 
 9088  T(2,|45)|4α=↓|43, T(3,|45)|4α=↓|44, T(4,|45)|4α=↓|43, 
 9091  |πso|'{A6}|εT|β5|4α=↓|4|f1|d35|)(1|4α+↓|42|4α+↓|43|4α+↓|44|4
 9092  α+↓|43)|4α=↓|42|f3|d35|).|;{A6}|π!|4|4|4|4In 
 9094  order to estimate |εT|βn |πfor large |εn, |πlet 
 9102  us _rst try an approximation suggested by R. 
 9110  W. Floyd: We might assume that, for |ε0|4|¬E|4k|4|¬W|4n, 
 9118  n |πis essentially ``random'' modulo |εk, |πso 
 9125  that we can set|'{A6}|εT|βn|4|¬V|41|4α+↓|4|(1|d2n|)|4(T|β0|4
 9129  α+↓|4T|β1|4α+↓|4|¬O|4|¬O|4|¬O|4α+↓|4T|βn|βα_↓|β1).|;
 9130  {A6}|πThen |εT|βn|4|¬V|4S|βn, |πwhere the sequence 
 9135  |ε|↔bS|βn|↔v |πis the solution to the recurrence 
 9142  relation|'{A6}|εS|β0|4α=↓|40,!!S|βn|4α=↓|41|4α+↓|4|(1|d2n|)|
 9143  4(S|β0|4α+↓|4S|β1|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4S|βn|βα_↓|β1
 9143  ),!!n|4|¬R|41.|J!(20)|;{A6}|π(This approximation 
 9146  is analogous to the ``lattice-point model'' used 
 9153  to investigate Algorithm B in Section 4.5.2.)|'
 9160  !|4|4|4|4The recurrence (20) is readily solved 
 9166  by the use of generating functions. A more direct 
 9175  way to solve it, analogous to ur solution of 
 9184  the lattice-point model, is by noting that|'{A6}|h|εS|βn|βα+
 9191  ↓|β1|4|∂α=↓|41|4α+↓|4n|4α+↓|41|4(n(S|βn|4α_↓|41)|4α+↓|4S|βn)
 9191  |4α=↓|4S|βn|4α+↓|4n|4α+↓|41|4;|E|n|;| S|βn|βα+↓|β1|4|Lα=↓|41
 9192  |4α+↓|4|(1|d2n|4α+↓|41|)|4(S|β0|4α+↓|4S|β1|4α+↓|1|1|¬O|4|¬O|
 9192  4|¬O|1|1α+↓|4S|βn|βα_↓|β1|4α+↓|4S|βn)>{A6}|Lα=↓|41|4α+↓|4|(1
 9193  |d2n|4α+↓|41|)|4{H12}({H10}n(S|βn|4α_↓|41)|4α+↓|4S|βn{H12})|
 9193  4{H10}α=↓|4S|βn|4α+↓|4|(1|d2n|4α+↓|41|)|4;>{A6}|πhence 
 9195  |εS|βn |πis 1|4α+↓|4|f1|d32|)|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4
 9197  1/|εn|4α=↓|4H|βn, |πa harmonic number. Therefore 
 9202  the approximation |εT|βn|4|¬V|4S|βn |πwould tell 
 9207  us that |εT|βn|4|¬V|4|πln|4|εn|4α+↓|4O(1).|'|π!|4|4|4|4Compa
 9210  rison of this approximation with tables of the 
 9218  true value of |εT|βn |πshow, however, that ln 
 9226  |εn |πis too large; |εT|βn |πdoes not grow this 
 9235  fast. One way to account for the fact that this 
 9245  approximation is too pessimistic is to observe 
 9252  that the average value of |εn|4|πmod|4|εk |πis 
 9259  less than the average value of |ε|f1|d32|)k, 
 9266  |πin the range |ε1|4|¬E|4k|4|¬E|4n:|'{A6}|(1|d2n|)|4|↔k|uc|)
 9270  1|¬Ek|¬En|)|4(n|4|πmod|4|εk)|4|∂α=↓|4|(1|d2n|)|4|↔k|)1|¬Eq|¬
 9270  En|d|"ln/(qα+↓1)|"L|¬Wk|¬E|"ln/q|"L|)|4(n|4α_↓|4qk)|4α=↓|4n|
 9270  4α_↓|4|(1|d2n|)|4|↔k|uc|)1|¬Eq|¬En|)|4q|↔a|↔a|(|"ln/q|"L|4α+
 9270  ↓|41|d52|)|↔s|4α_↓|4|↔a|(|"ln/(q|4α+↓|41)|"L|4α+↓|41|d52|)|↔
 9270  s|↔s|;{A6}|Lα=↓|4n|4α_↓|4|(1|d2n|)|4|↔k|uc|)1|¬Eq|¬En|)|4|↔a
 9271  |(|"ln/q|"L|4α+↓|41|d52|)|↔s|4α=↓|4|↔a1|4α_↓|4|(|≤p|g2|d212|
 9271  )|↔sn|4α+↓|4O(|πlog|4|εn)|J!(21)>{A6}{H12}({H10}cf. 
 9273  exercise 4.5.2<10(c){H12}){H10}. This is only 
 9278  about .1775|εn, |πnot |ε.25n, |πso the value 
 9285  of |εn |πmod|4|εk |πtends to be smaller than 
 9293  the above model predicts; Euclid's algorithm 
 9299  works fsster than we might expect.|'|H*?*?{U0}{H10L12M29}58320
folio 448 galley 13
 9305  #Computer Programming!(A.-W./Knuth)!f. 448!ch. 
 9308  4!g. 13|'{A12}|≡A|≡c|≡o|≡n|≡t|≡i|≡n|≡u|≡o|≡u|≡s 
 9311  |≡m|≡o|≡d|≡e|≡l|≡.|9|4The behavior of Euclid's 
 9315  algorithm with |εv|4α=↓|4N |πis essentially determined 
 9321  by the behavior of the regular continued fraction 
 9329  process when |εX|4α=↓|40/N, 1/N,|4.|4.|4.|4,|4(N|4α_↓|41)/N.
 9332   |πAssuming that |εN |πis very large, we are 
 9341  led naturally to a study of regular continued 
 9349  fractions when |εX |πis a real number uniformly 
 9357  distributed in [0,|41). Therefore let us de_ne 
 9364  the distribution function|'{A6}|εF|βn(x)|4α=↓|4|πprobability
 9367  |4that|4|εX|βn|4|¬E|4x,!!|πfor!!|ε0|4|¬E|4x|4|¬E|41,|J!(22)|
 9367  ;{A6}|πgiven a uniform distribution of |εX|4α=↓|4X|β0. 
 9374  |πBy the de_nition of regular continued fractions, 
 9381  we have |εF|β0(x)|4α=↓|4x, |πand|'{A6}|h|εF|βn|βα+↓|β1(x)|4|
 9385  Lα=↓|4!!(|πprobability|4that|4|ε1/(k|4α+↓|4x)|4α+↓|4X|βn|4α+
 9385  ↓|41/k)|E|n|;| F|βn|βα+↓|β1(x)|4|Lα=↓|4|↔k|uc|)k|¬R1|)|4(|πp
 9386  robability|4that|4|εk|4|¬E|41/X|βn|4|¬E|4k|4α+↓|4x)>
 9387  {A6}|Lα=↓|4|↔k|uc|)k|¬R1|)|4(|πprobability|4that|4|ε1/(k|4α+
 9387  ↓|4x)|4|¬E|4X|βn|4|¬E|41/k)>{A6}|Lα=↓|4|↔k|uc|)k|¬R1|)|4(F|β
 9388  n(1/k)|4α_↓|4F|βn(1/(k|4α+↓|4){H12})){H10}.|J!(23)>
 9389  {A6}|πIf the distributions |εF|β0(x),|4F|β1(x),|4.|4.|4. 
 9393  |πde_ned by these formulas approach a limiting 
 9400  distribution |εF|β|¬X(x)|4α=↓|4F(x), |πwe will 
 9404  have|'{A6}|εF(x)|4α=↓|4|↔k|uc|)k|¬R1|)|4{H12}({H10}F(1/k)|4α
 9405  _↓|4F{H12}({H10}1/(k|4α+↓|4x){H12})){H10}.|J!(24)|;
 9406  {A6}|π!|4|4|4|4One function which satis_es this 
 9411  relation is |εF(x)|4α=↓|4|πlog|ε|βb|4(1|4α+↓|4x), 
 9414  |πfor any base |εb|4|¬Q|41; |πsee exercise 19. 
 9421  The further condition |εF(1)|4α=↓|41 |πimplies 
 9426  that we should take |εb|4α=↓|42. |πThus it is 
 9434  reasonable to make a guess that |εF(x)|4α=↓|4|πlg|β2|4(1|4α+
 9440  ↓|4|εx), |πand that |εF|βn(x) |πapproaches this 
 9446  behavior.|'!|4|4|4|4We might conjecture, for 
 9451  example, that |εF(|f1|d32|))|4α=↓|4|πlg(|f3|d32|))|4|¬V|40.5
 9453  8496; let us see how close |εF|βn(|f1|d32|)) 
 9460  |πcomes to this value for small |εn. |πWe have 
 9469  |εF|β0(|f1|d32|)), |πand|'|ε|h|εF|β2(|9)|4|∂α=↓|4!!m|9|1|12m
 9471  |4α+↓|42|4α_↓|43m|4α_↓|42|4α+↓|44m|4α+↓|42|4α+↓|45m|4α+↓|42|
 9471  4α+↓|4.|4.|4.|9|1|1|E|n|;|ε| F|β1(|f1|d32|))|4|Lα=↓|4|(1|d21
 9472  |)|4α_↓|4|(1|d21|4α+↓|4|f1|d32|)|)|4α+↓|4|(1|d22|)|4α_↓|4|(1
 9472  |d22|4α+↓|4|f1|d32|)|)|4α+↓|1|1|¬O|4|¬O|4|¬O>
 9473  {A5}|ε|Lα=↓|42|↔a|(1|d22|)|4α_↓|4|(1|d23|)|4α+↓|4|(1|d24|)|4
 9473  α_↓|4|(1|d25|)|4α+↓|1|1|¬O|4|¬O|4|¬O|4|↔s|4α=↓|42(1|4α_↓|4|π
 9473  ln|42)|4|¬V|40.6137;>|ε{A5}| |εF|β2(|f1|d32|))|4|Lα=↓|4|↔k|u
 9474  c|)m|¬R1|)|4|(2|d2m|)|4|↔a|(1|d22m|4α+↓|42|)|4α_↓|4|(1|d23m|
 9474  4α+↓|42|)|4α+↓|4|(1|d24m|4α+↓|42|)|4α_↓|4|(1|d25m|4α+↓|42|)|
 9474  4α+↓|1|1|¬O|4|¬O|4|¬O|1|↔s>|ε{A4}|Lα=↓|4|↔k|uc|)m|¬R1|)|4|(2
 9475  |d2m|g2|)|4|↔a|(1|d22|)|4α_↓|4|(1|d23|)|4α+↓|4|(1|d24|)|4α_↓
 9475  |1|1|¬O|4|¬O|4|¬O|↔s>|ε{A4}|L|4|4|4|4|4α_↓|4|↔k|uc|)m|¬R1|)|
 9476  4|(4|d2m|)|4|↔a|(1|d22m(2m|4α+↓|42)|)|4α_↓|4|(1|d23m(3m|4α+↓
 9476  |42)|4α+↓|1|1|¬O|4|¬O|4|¬O|↔s>{A5}|ε|Lα=↓|4|f1|d23|)|≤p|g2(1
 9477  |4α_↓|4|πln|42)|4α_↓|4|↔k|uc|)m|¬R1|)|4|(4S|βm|d2m|g2|)|4,>
 9478  {A6}|πfor some positive constant |εA. |πThe error 
 9485  term was improved to |εO(e|gα_↓|gA|gn) |πby Paul 
 9492  L|=1evy shortly afterward [|εBull. Soc. Math. 
 9498  de France |≡5|≡7 (1929), 178<194];|≤⊂ |πbut Gauss's 
 9505  problem, namely to _nd the asymptotic behavior 
 9512  of |εF|βn(x)|4α_↓|4|πlg(1|4α+↓|4|εx), |πwas not 
 9516  really resolved until 1974, when Eduard Wirsing 
 9523  published a beautiful analysis of the situation 
 9530  [|εActa Arithmetica |≡2|≡4 (1974), 507<528]. 
 9535  |πWe shall study the simplest aspects of Wirsing's 
 9543  approach here, since his method appears to be 
 9551  the best way to handle problems of this type.|'
 9560  !|4|4|4|4If |εG |πis any function of |εx |πde_ned 
 9568  for |ε0|¬E|4x|4|¬E|41, |πlet |εSG |πbe the function 
 9575  de_ned by|'{A6}|εSG(x)|4α=↓|4|↔k|uc|)k|¬R1|)|4|↔aG|↔a|(1|d2k
 9577  |)|↔s|4α_↓|4G|↔a|(1|d2k|4α+↓|4x|)|↔s|↔s.|J!(25)|;
 9578  {A6}|πThus, |εS |πis an operator which changes 
 9585  one function into another. In particular, by 
 9592  (23) we have |εF|βn|βα+↓|β1(x)|4α=↓|4SF|βn(x), 
 9596  |πhence|'{A6}|εF|βn|4α=↓|4S|gnF|β0.|J!(26)|;|π{A6}(In 
 9599  this discussion |εF|βn |πstands for a distribution 
 9606  function, |εnot |πfor a Fibonacci number.) Note 
 9613  that |εS |πis a ``linear operator''; i.e., |εS(cG)|4α=↓|4c(S
 9620  G) |πfor constants |εc, |πand |εS(G|β1|4α+↓|4G|β2)|4α=↓|4SG|
 9625  β1|4α+↓|4SG|β2.|'!|4|4|4|4|πNow if |εG |πhas 
 9630  a bounded _rst derivative, we can di=erentiate 
 9637  (25) term by term to show that|'{A6}|ε(SG)|¬S(x)|4α=↓|4|↔k|u
 9644  c|)k|¬R1|)|4|(1|d2(k|4α+↓|4x)|g2|)|4G|¬S|↔a|(1|d2k|4α+↓|4x|)
 9644  |↔s,|J!(27)|;{A6}#####|'{A3}{H8}|≤⊂|4|πAn exposition 
 9648  of L|=1evy's interesting proof appeared in the 
 9655  _rst edition of this book.|'{A6}{H10}hence |εSG 
 9662  |πalso has a bounded _rst derivative. (Term-by-term 
 9669  di=erentiation of a convergent series is justi_ed 
 9676  when the series of derivatives is uniformly convergent; 
 9684  cf. K. Knopp, |εTheory and Application of In⊂nite 
 9692  series (|πGlasgow: Blackie, 1951), |9|147.)|'
 9697  !|4|4|4|4Let |εH|4α=↓|4SG, |πand let |εg(x)|4α=↓|4(1|4α+↓|4x
 9701  )G|¬S(x), h(x)|4α=↓|4(1|4α+↓|4x)H|¬S(x). |πIt 
 9704  follows that|'{A6}|εh(x)|4α=↓|4|↔k|uc|)k|¬R1|)|4|(1|4α+↓|4x|
 9706  d2(k|4α+↓|4x)|g2|)|4|↔a1|4α+↓|4|(1|d2k|4α+↓|4x|↔s|gα_↓|g1g|↔
 9706  a|(1|d2k|4α+↓|4x|↔s|4α=↓|↔k|uc|)k|¬R1|)|1|1|↔a|(k|d2k|4α+↓|4
 9706  1|4α+↓|4x)|4α_↓|4|(k|4α_↓|41|d2k|4α+↓|4x|)|↔sg|↔a|(1|d2k|4α+
 9706  ↓|4x|↔s|)|↔s.|;{A6}|πIn other words, |εh|4α=↓|4Tg, 
 9711  |πwhere |εT |πis the linear operator de_ned by|'
 9719  {A6}|εTg(x)|4α=↓|4|↔k|uc|)k|¬R1|)|4|↔a|(k|d2k|4α+↓|41|4α+↓|4
 9719  x|)|4α_↓|4|(k|4α_↓|41|d2k|4α+↓|4x|)|↔sg|↔a|(1|d2k|4α+↓|4x|)|
 9719  ↔s.|J!(28)|;!|4|4|4|4Continuing, we see that 
 9724  if |εg |πhas a bounded _rst derivative, we can 
 9733  di=erentiate term by term to show that |εTg |πdoes 
 9742  also:|'{A6}(Tg)|¬S(x)|4|∂α=↓|4α_↓|4|↔k|uc|)k|¬R1|)|4|↔a|(k|d
 9743  2(k|4α+↓|41|4α+↓|4x)|g2|)|4|↔ag|↔a|(1|d2k|4α+↓|4x|)|↔s|4α_↓|
 9743  4g|↔a|(1|d2k|4α+↓|41|4α+↓|4x|)|↔s|↔s|4α+↓|4|(1|4α+↓|4x|d2(k|
 9743  4α+↓|4x)|g3(k|4α+↓|41|4α+↓|4x)|)|4g|¬S|↔a|(1|d2k|4α+↓|4x|)|↔
 9743  s|↔s.|;{A6}|Lα=↓|4α_↓|4|↔k|uc|)k|¬R1|)|4|↔a|(k|d2(k|4α+↓|41|
 9744  4α+↓|4x)|g2|)|4|↔ag|↔a|(1|d2k|4α+↓|4x|)|↔s|4α_↓|4g|↔a|(1|d2k
 9744  |4α+↓|41|4α+↓|4x|)|↔s|↔s|4α+↓|4|(1|4α+↓|4x|d2(k|4α+↓|4x)|g3(
 9744  k|4α+↓|41|4α+↓|4x)|)|4g|¬S|↔a|(1|d2k|4α+↓|4x|)|↔s|↔s.>
 9745  {A6}|πThere is consequently a third linear operator, 
 9752  |εU, |πsuch that |ε(Tg)|¬S|4α=↓|4α_↓U(g)|¬S), 
 9756  |πnamely|'{A6}|εU|≤'(x)|4α=↓|4|↔k|uc|)k|¬R1|)|4|↔a|(k|d2(k|4
 9757  α+↓|41|4α+↓|4x)|g2|)|4|↔j|ur1/(kα+↓x)|)1/(kα+↓1α+↓x)|)|4|≤'(
 9757  t)dt|4α+↓|4|(1|4α+↓|4x|d2(k|4α+↓|4x)|g3(k|4α+↓|41|4α+↓|4x)|)
 9757  |4|≤'|↔a|(1|d2k|4α+↓|4x|)|↔s|↔s.|;(29)|?{A6}|H*?*?*?*?{U0}{H10L1
folio 451 galley 14
 9759  2M29}58320#Computer Programming!(A.-W./Knuth)!ch. 
 9761  4!f. 451!g. 14|'{A12}|πWhat is the relevance 
 9768  of all this to our problem? Well, if we set|'
 9778  {A6}|h|εf|βn(x)|4|∂α=↓|4(1|4α+↓|4x)F|βn(x)|4α=↓|4|πln|42|4(1
 9778  |4α+↓|4|εR|βn(|πlg(1|4α+↓|4|εx))),|E|n|;| F|βn(x)|4|Lα=↓|4|π
 9779  lg|ε(1|4α+↓|4x)|4α+↓|4R|βn{H12}({H10}|πlg(1|4α+↓|4x){H12}){H
 9779  10},|J!(30)>{A5}| |εf|βn(x)|4|Lα=↓|4(1|4α+↓|4x)F|1|ur|↔0|)n|
 9780  )(x)|4α=↓|4|(1|d2|πln|42|)|4|ε{H12}({H10}1|4α+↓|4R|ur|↔0|)n|
 9780  ){H12}({H10}|πlg(1|4α+↓|4x){H12})){H10},|J!(31)>
 9781  {A6}we have|'{A6}|ε| f|1|ur|↔0|)n|)(x)|4|Lα=↓|4R|βn(|πlg(1|4
 9783  α+↓|4|εx))/(|πln|42)|g2(1|4α+↓|4|εx);|J!(32)>
 9784  {A6}|πthe e=ect of the lg(1|4α+↓|4|εx) |πterm 
 9790  disappears, after these transformations. Furthermore 
 9795  since |εF|βn|4α=↓|4S|gnF|β0 |πwe have |εf|βn|4α=↓|4T|gnf|β0 
 9800  |πand |εf|1|ur|↔0|)n|)|4α=↓|4(α_↓1)|gnU|gnf|1|ur|↔0|)0|), 
 9802  |πand |εF|βn |πand |εf|βn |πhave bounded derivatives, 
 9809  by induction on |εn. |πThus (32) becomes|'{A6}|ε(|→α_↓1)|gnR
 9816  |βn{H12}(|π{H10}lg(1|4α+↓|4|εx){H12})|4{H10}α=↓|4(1|4α+↓|4x)
 9816  |πln|42)|g2|εU|gnf|1|ur|↔0|)0|)(x).|J!(33)|;{A6}|π{H10L12M29
 9817  }Now |εF|β0(x)|4α=↓|4x, f|β0(x)|4α=↓|41|4α+↓|4x, 
 9820  |πand |εf|β|1|ur|↔0|)0|)(x) |πis the constant 
 9825  function 1. We are going to show that the operator 
 9835  |εU|gn |πtakes the constant function into a function 
 9843  with very small values, hence |ε|¬GR|βn(x)|¬G 
 9849  |πmust be very small for |ε0|4|¬E|4x|4|¬E|41. 
 9855  |πFinally we can clinch the argument by showing 
 9863  that |εR|βn(x) |πitself is small: Since |εR|βn(0)|4α=↓|4R|βn
 9869  (1)|4α=↓|41, |πit follows from a well-known interpolation 
 9876  formula (cf. exercise 4.6.4<15 with |εx|β0|4α=↓|40, 
 9882  x|β1|4α=↓|4x, x|β2|4α=↓|41) |πthat|'{A6}|εR|βn(x)|4α=↓|4α_↓|
 9885  (x(1|4α_↓|4x)|d22|)|4R|=|βn|g|¬C(|≤j(x))|J!(34)|;
 9886  |πfor some function |ε|≤j(x), |πwhere |ε0|4|¬E|4|≤j(x)|4|¬E|
 9891  41 |πwhen |ε0|4|¬E|4x|4|¬E|41.|'|π!|4|4|4|4Thus 
 9895  everything hinges on our being able to prove 
 9903  that |εU|gn |πproduces small function values, 
 9909  where |εU |πis the linear operator de_ned in 
 9917  (29). Note that |εU |πis a |εpositive |πoperator 
 9925  in the sense that |εU|≤'(x)|4|¬R|40 |πfor all 
 9932  |εx |πif |ε|≤'(x)|4|¬R|40 |πfor all |εx. |πIt 
 9939  follows that |εU |πis order-preserving: If |ε|≤'|β1(x)|4|¬E|
 9945  4|≤'|β2(x) |πfor all |εx |πthen |εU|≤'|β1(x)|4|¬E|4U|≤'|β2(x
 9950  ) |πfor all |εx.|'|π!|4|4|4|4One way to exploit 
 9958  this property is to _nd a function |ε|≤' |πfor 
 9967  which we can calculate |εU|≤' |πexactly and to 
 9975  use constant multiples of this function to bound 
 9983  the ones were really interested in. First let's 
 9991  look for a function |εg |πsuch that |εTg |πis 
10000  easy to compute. If we consider functions de_ned 
10008  for all |εx|4|¬R|40 |πinstead of only on [0,|41], 
10016  it is easy to remove the summation from (25) 
10025  by observing that|'{A6}|εSG(x|4α+↓|41)|4α_↓|4SG(x)|4α=↓|4G|↔
10028  a|(1|d21|4α+↓|4x|)|↔s|4α_↓|4|uw|πlim|uc|)|εk|¬M|¬X|)|4G|↔a|(
10028  1|d2k|4α+↓|4x|)|↔s|4α=↓|4G|↔a|(1|d21|4α+↓|4x|)|↔s|4α_↓|4G(0)
10028  |;|da3}(35)|?{A6}|πwhen |εG |πis continuous. 
10034  Since |εT{H12}({H10}(1|4α+↓|4x)G|¬S{H12}){H10}|4α=↓|4(1|4α+↓
10035  |4x)(SG)|¬S, |πit follows (see exercise 20) that|'
10042  {A6}|ε|(Tg(x)|d2x|4α+↓|41|)|4α_↓|4|(Tg(x|4α+↓|41)|d2x|4α+↓|4
10042  2|)|4α=↓|4|↔a|(1|d2x|4α+↓|41|)|4α_↓|4|(1|d2x|4α+↓|42|)|↔sg|↔
10042  a|(1|d21|4α+↓|4x|)|↔s.|J!(36)|;{A6}|πIf we set 
10046  |εTg(x)|4α=↓|41/(x|4α+↓|41), |πwe _nd that the 
10051  corresponding value of |εg(x) |πis |ε1|4α+↓|4x|4α_↓|41/(1|4α
10056  +↓|4x). |πLet |ε|≤'(x)|4α=↓|4g|¬S(x)|4α=↓|41|4α+↓|41/(1|4α+↓
10058  |4x)|g2, |πso that |εU|≤'(x)|4α=↓|4α_↓(Tg)|¬S(x)|4α=↓|41/(1|
10061  4α+↓|4x)|g2; |πthis is the function |ε|≤' |πwe 
10068  have been looking for.|'!|4|4|4|4For this choice 
10075  of |ε|≤' |πwe have |ε2|4|¬E|4|≤'(x)/U|≤'(x)|4α=↓|4(1|4α+↓|4x
10079  )|g2|4α+↓|41|4|¬E|45 |πfor |ε0|4|¬E|4x|4|¬E|41, 
10082  |πhence|'{A6}|ε|f1|d35|)|≤'|4|¬E|4U|≤'|4|¬E|4|f1|d32|)|≤'.|;
10084  {A6}|πBy the positivity of |εU |πand |ε|≤' |πwe 
10092  can apply |εU |πto this inequality again, obtaining 
10100  |ε|f1|d325|)|≤'|4|¬E|4|f1|d35|)U|≤'|4|¬E|4U|g2|≤'|4|¬E|4U|≤'
10100  |4|≤E|4|f1|d34|)|≤'; |πand after |εn|4α_↓|41 
10104  |πapplications we have|'{A6}|ε5|gα_↓|gn|≤'|4|¬E|4U|gn|≤'|4|¬
10107  E|42|gα_↓|gn|≤'|J!(37)|;{A6}|πfor this particular 
10111  |ε|≤'. |πLet |ε|≤c(x)|4α=↓|4f|1|ur|↔0|)0|)(x)|4α=↓|41 
10114  |πbe the constant function; then for |ε0|4|¬E|4x|4|¬E|41 
10121  |πwe have |ε|f5|d34|)|≤c|4|¬E|4|≤'|4|¬E|42|≤c, 
10124  |πhence|'{A6}|ε|f5|d38|)5|gα_↓|gn|≤c|4|¬E|4|f1|d32|)5|gα_↓|g
10125  n|≤'|4|¬E|4|f1|d32|)U|gn|≤'|4|¬E|4U|gn|≤c|4|¬E|4|f4|d35|)U|g
10125  n|≤'|4|¬E|4|f4|d35|)2|gα_↓|gn|≤'|4|¬E|4|f8|d35|)2|gα_↓|gn|≤c
10125  .|;{A6}It follows by (33) that|'{A6}|ε|f5|d38|)(|πln|4|ε2)|g
10131  25|gα_↓|gn|4|¬E|4(|→α_↓1)|gnR|βn(x)|4|¬E|4|f16|d35|)(|πln|4|
10131  ε2)|g22|gα_↓|gn,!0|4|¬E|4x|4|¬E|41,|;{A6}|πhence 
10133  by (30) and (34) we have proved the following 
10142  result:|'{A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡W|≡.|9|4|εF|βn(x)|4α=↓
10144  |4|πlg(1|4α+↓|4|εx)|4α+↓|40(2|gα_↓|gn). In fact, 
10147  F|βn(x)|4α_↓|4|πlg(1|4α+↓|4|εx) lies between 
10150  |f5|d316|)(|→α_↓1)|gn|gα+↓|g15|gα_↓|gn{H12}({H10}|πln(1|4α+↓
10150  |4|εx){H12}){H10}(|πln|42/1|4α+↓|4|εx) and |f8|d35|)(|→α_↓1)
10152  |gn|gα+↓|g12|gα_↓|gn{H12}({H10}|πln(1|4α+↓|4|εx){H12}){H10}(
10152  |πln|42/1|4α+↓|4|εx), for 0|4|¬E|4x|4|¬E|41.|'
10155  {A12}|πWith a slightly di=erent choice of |ε|≤', 
10162  |πwe can obtain tighter bounds. (see exercise 
10169  21). In fact, Wirsing went much further in his 
10178  paper, proving that|'{A6}|εF|βn(x)|4α=↓|4|πlg(1|4α+↓|4|εx)|4
10181  α+↓|4(|→α_↓|≤l)|gn|≤<(x)|4α+↓|4O(x(1|4α_↓|4x)(|≤l|4α_↓|40.03
10181  1)|gn),|J!(38)|;{A6}|πwhere|'{A6}|h|ε|≤l|4|∂α=↓|4|"C3,|43,|4
10183  2,|42,|43,|413,|41,|4174,|41,|41,|41,|42,|42,|42,|41,|41,|41
10183  ,|42,|42,|41,|4.|4.|4.|"C|E|n|;| |≤l|4|Lα=↓|40.30366|930028|
10184  998732|965860|4.|4.|4.>{A4}|Lα=↓|4|"C3,|43,|42,|42,|43,|413,
10185  |41,|4174,|41,|41,|41,|42,|42,|42,|41,|41,|41,|42,|42,|41,|4
10185  .|4.|4.|"C|J!(39)|;{A6}|πis an apparently new 
10190  fundamental constant, and where |ε|≤< |πis an 
10197  interesting function which is analytic in the 
10204  entire complex plane except for the negative 
10211  real axis from |→α_↓1 to |→α_↓|¬X. Wirsing's 
10218  function satis_es |ε|≤<(0)|4α=↓|4|≤<(1)|4α=↓|40, 
10221  |≤<|¬S(0)|4|¬W|40, |πand |εS|≤<|4α=↓|4|→α_↓|≤l|≤<; 
10224  |πthus by (35) it satis_es the identity|'{A6}|ε|≤<(z)|4α_↓|4
10231  |≤<(z|4α+↓|41)|4α=↓|4|(1|d2|≤l|)|4|≤<|↔a|(1|d21|4α+↓|4z|)|↔s
10231  .|J!(40)|;{A6}|πFurthermore, Wirsing demonstrated 
10235  that|'{A6}|ε|≤<|↔aα_↓|(u|d2v|)|4α+↓|4|(i|d2N|)|↔s|4α=↓|4c|≤l
10236  |gα_↓|gn|4|πlog|4|εN|4α+↓|4O(1)!|πas!|εN|4|¬M|4|¬X,|J!(41)|;
10237  {A6}|πwhere |εc |πis a constant and |εn|4α=↓|4T(u,|4v) 
10244  |πis the number of iterations when Euclid's algorithm 
10252  is applied to the integers |εu|4|¬Q|4v|4|¬Q|40.|'
folio 460 galley 15
11314  2M29}58320#Computer Programming!(A.-W./Knuth)!ch. 
11316  4!f. 460!g. 15|'{A12}|π|≡F|≡r|≡o|≡m |≡c|≡o|≡n|≡t|≡i|≡n|≡u|≡o
11320  |≡u|≡s |≡t|≡o |≡d|≡i|≡s|≡c|≡r|≡e|≡t|≡e|≡.|9|4We 
11323  have now derived results about the probability 
11330  distributions for continued fractions when |εX 
11336  |πis a real number uniformly distributed in the 
11344  interval [0,|41). But since a real number is 
11352  rational with probability zero (almost all numbers 
11359  are irrational), these results do not apply directly 
11367  to Euclid's algorithm. Before we can apply Theorem 
11375  W to our problem, some technicalities must be 
11383  overcome. Consider the following observation 
11388  based on elementary measure theory:|'{A12}|≡L|≡e|≡m|≡m|≡a 
11394  |≡M|≡.|9|4|εLet I|β1, I|β2,|4.|4.|4.|4,|4J|β1,|4J|β2,|4.|4.|
11396  4. be pairwise disjoint intervals contained in 
11403  the interval [0,|41), and let |λI|4α=↓|4|↔e|βk|β|¬R|β1|4I|βk
11408  , |λJ|4α=↓|4|↔e|βk|β|¬R|β1|4J|βk. Assume that 
11412  |λK|4α=↓|4[0,|41]|4α_↓|4(|λI|4|↔q|4|λJ) has measure 
11415  zero. Let P|βn be the set |¬T0/n, 1/n,|4.|4.|4.|4,|4(n|4α_↓|
11422  41)/n|¬Y. Then|'{A6}|uw|πlim|uc|)|εn|¬M|¬X|)|4|(|¬F|λI|4|↔Q|
11424  4P|βn|¬F|d2n|)|4α=↓|4|≤m(|λI).|J!(42)|;{A6}|πHere 
11426  |ε|≤m(|λI) |πis the Lebesgue measure of |ε|λI, 
11433  |πnamely, |ε|¬K|βk|β|¬R|β1 |πlength |ε(I|βk); 
11437  |πand |ε|¬F|λi|4|↔Q|4P|βn|¬F |πdenotes the number 
11442  of elements of |ε|λI|4|↔Q|4P|βn.)|'{A12}Proof.|9|4|πLet 
11447  |ε|λI|βN|4α=↓|4|↔e|ur|)1|¬Ek|¬EN|)|4I|βk, |λJ|βN|4α=↓|4|↔e|u
11448  r|)1|¬Ek|¬EN|)|4J|βk. |πGiven |ε|≤e|4|¬Q|40, 
11451  |π_nd |εN |πlarge enough so that |ε|≤m(|λI|βN)|4α+↓|4|≤m(|λJ
11457  |βN)|4|¬R|41|4α_↓|4|≤e, |πand let|'{A6}|ε|λK|βN|4α=↓|4|λK|4|
11460  ↔q|4|↔e|βk|β|¬Q|βN|4I|βk|4|↔q|4|↔e|βk|β|¬Q|βN|4J|βk.|;
11461  {A6}|πIf |εI |πis an interval, having any of 
11469  the forms |ε(a,|4b) |πor |ε[a,|4b) |πor |ε(a,|4b] 
11476  |πor |ε[a,|4b], |πit is clear that |ε|≤m(I)|4α=↓|4b|4α_↓|4a 
11483  |πand|'{A6}|εn|≤m(I)|4α_↓|41|4|¬E|4|¬FI|4|↔Q|4P|βn|¬F|4|¬E|4
11484  n|≤m(I)|4α+↓|41.|;{A6}|πNow let |εr|βn|4α=↓|4|¬F|λI|βN|4|↔Q|
11487  4P|βn|¬F, s|βn|4α=↓|4|¬F|λJ|βN|4|↔Q|4P|βn|¬F, 
11489  t|βn|4α=↓|4|¬F|λK|βN|4|↔Q|4P|βn|¬F; |πwe have|'
11492  {A6}|εr|βn|4α+↓|4s|βn|4α+↓|4t|βn|4α=↓|4n;|;{A6}n|≤m(|λI|βN)|
11493  4α_↓|4N|4|¬E|4r|βn|4|¬E|4n|≤m(|λI|βN)|4α+↓|4N;|;
11494  {A3}n|≤m(|λI|βN)|4α_↓|4N|4|¬E|4s|βn|4|¬E|4n|≤m(|λJ|βN)|4α+↓|
11494  4N.|;{A6}|πHence|'{A6}|ε|≤m(|λI)|4α_↓|4|(N|d2n|)|4α_↓|4|≤e|4
11496  |¬E|4|≤m(|λI|βN)|4α_↓|4|(N|d2n|)|4|¬E|4|(r|βn|d2n|)|4|¬E|4|(
11496  r|βn|4α+↓|4t|βn|d2n|)|'{A4}α=↓|41|4α_↓|4|(s|βn|d2n|)|4|¬E|41
11497  |4α_↓|4|≤m(|λJ|βN)|4α+↓|4|(N|d2n|)|4|¬E|4|≤m(|λI)|4α+↓|4|(N|
11497  d2n|)|4α+↓|4|≤e.|?{A6}|πThis holds for all |εn 
11503  |πand for all |ε|≤e; |πhence lim|ε|βn|β|¬M|β|¬X|4r|βn/n|4α=↓
11508  |4|≤m(|λI).|'{A12}|π!|4|4|4|4Exercise 25 shows 
11512  that Lemma M is not trivial, in the sense that 
11522  some rather restrictive hypotheses are needed 
11528  to prove (42).|'{A12}|≡D|≡i|≡s|≡t|≡r|≡i|≡b|≡u|≡t|≡i|≡o|≡n 
11532  |≡o|≡f |≡p|≡a|≡r|≡t|≡i|≡a|≡l |≡q|≡u|≡o|≡t|≡i|≡e|≡n|≡t|≡s|≡.|
11534  9|4Now we can put Theorem W and Lemma M together 
11544  to derive some solid facts about Euclid's algorithm.|'
11552  {A12}|≡T|≡h|≡e|≡o|≡r|≡e|≡m |≡E|≡.|9|4|εLet n 
11555  and k be positive integers, and let p|βk(a,|4n) 
11563  be the probability that the (k|4α+↓|41)st quotient 
11570  A|βk|βα+↓|β1 in Euclid's algorithm is equal to 
11577  a, when v|4α=↓|4n and u is chosen at random. 
11586  Then|'{A6}|uw|πlim|uc|)|εn|¬M|¬X|)|4p|βk(a,|4n)|4α=↓|4F|βk|↔
11587  a|(1|d2a|)|↔s|4α_↓|4F|βk|↔a|(1|d2a|4α+↓|41|)|↔s,|;
11588  {A6}where F|βk(x) is the distribution function 
11594  (21).|'{A12}Proof.|9|4|πThe set |ε|λI |πof all 
11600  |εX |πin [0,|41) for which |εA|βk|βα+↓|β1|4α=↓|4a 
11606  |πis a union of disjoint intervals, and so is 
11615  the set |ε|λJ |πof all |εX |πfor which |εA|βk|βα+↓|β1|4|=|↔6
11623  α=↓|4a. |πLemma M therefore applies, with |ε|λK 
11630  |πthe set of all |εX |πfor which |εA|βk|βα+↓|β1 
11638  |πis unde_ned. Furthermore, |εF|βk(1/a)|4α_↓|4F|βk(1/(a|4α+↓
11641  |41){H12}) {H10}|πis the probability that |ε1/(a|4α+↓|41)|4|
11646  ¬W|4X|βk|4|¬E|41/a, |πwhich is |ε|≤m(|λI), |πthe 
11651  probability that |εA|βk|βα+↓|β1|4α=↓|4a.|'{A12}|π!|4|4|4|4As
11654   a consequence of Theorems E and W, we can say 
11665  that a quotient equal to |εa |πoccurs with the 
11674  approximate probability|'{A6}|πlg(1|4α+↓|41/|εa)|4α_↓|4|πlg{
11676  H12}({H10}1|4α+↓|41/(|εa|4α+↓|41){H12})|4{H10}α=↓|4|πlg{H12}
11676  ({H10}(|εa|4α+↓|41)|g2/({H10}(a|4α+↓|41)|g2|4α_↓|41){H12}){H
11676  10}.|;{A6}|πThus|'{A6}|∂|h|πa|4|1quotient|4|1of|4|11|4|1occu
11678  rs|4|1about|4|1lg(|f23|d323|))|4|∂α=↓|433.333|4|1percent|4|1
11678  of|4|1the|4|1time;|E|n|;|La|4|1quotient|4|1of|4|11|4|1occurs
11679  |4|1about|4|1lg(|f4|d33|))|Lα=↓|441.504|4|1percent|4|1of|4|1
11679  the|4|1time;>{A3}|La|4|1quotient|4|1of|4|12|4|1occurs|4|1abo
11680  ut|4|1lg(|f9|d38|))|Lα=↓|416.992|4|1percent|4|1of|4|1the|4|1
11680  time;>{A3}|La|4|1quotient|4|1of|4|13|4|1occurs|4|1about|4|1l
11681  g(|f16|d315|))|Lα=↓|4|9|19.311|4|1percent|4|1of|4|1the|4|1ti
11681  me;>{A3}|La|4|1quotient|4|1of|4|14|4|1occurs|4|1about|4|1lg(
11682  |f25|d324|))|Lα=↓|4|9|15.890|4|1percent|4|1of|4|1the|4|1time
11682  .>{A6}Actually, if Euclid's algorithm produces 
11688  the quotients |εA|β1,|4A|β2,|4.|4.|4.|4,|4A|βt, 
11691  |πthe nature of the proofs above will guarantee 
11699  this behavior only for |εA|βk |πwhen |εk |πis 
11707  comparatively small with respect to |εt; |πthe 
11714  values |εA|βt|βα_↓|β1, A|βt|βα_↓|β2,|4.|4.|4. 
11717  |πare not covered by this proof. But we can in 
11727  fact show that the distribution of the last quotients 
11736  |εA|βt|βα_↓|β1, A|βt|βα_↓|β2,|4.|4.|4. |πis essentially 
11740  the same as the _rst.|'!|4|4|4|4For example, 
11747  consider the regular continued fraction expansions 
11753  for the fractions whose denominator is 29:|'|Hβ*?*?*!*!*!*!*!*!*!*!*!*!*!
11760  
folio 461 galley 16
fol58  |H*?*?*?{U0}{H10L12M29}58320#Computer Programming!(A.-W./Knuth)
10259  !ch. 4!f. 461!g. 16|'{A12}|∂!!!!!!!|1|1|∂!!!!!!!!!|4|∂!!!!!!
10263  !!!!|4|1|∂!!!!!!!|1|1|1|∂|E|'|>|f1|d329|)|4α=↓|4|"C29|"C|'
10266  |f8|d329|)|4α=↓|4|"C3,|41,|41,|41,|42|"C|'|f15|d329|)|4α=↓|4
10267  |"C1,|41,|414|"C|'|f22|d329|)|4α=↓|4|"C1,|43,|47|"C|'
10269  >|>|f2|d329|)|4α=↓|4|"C14,|42|"C|'|f9|d329|)|4α=↓|4|"C3,|44,
10272  |42|"C|'|f16|d329|)|4α=↓|4|"C1,|41,|44,|43|"C|'
10274  |f23|d329|)|4α=↓|4|"C1,|43,|41,|45|"C|'>|>|f3|d329|)|4α=↓|4|
10277  "C9,|41,|42|"c|'|f10|d329|)|4α=↓|4|"C2,|41,|49|"c|'
10279  |f17|d329|)|4α=↓|4|"C1,|41,|42,|42,|42|"C|'|f24|d329|)|4α=↓|
10280  4|"C1,|44,|41,|44|"C|'>|>|f4|d329|)|4α=↓|4|"C7,|44|"C|'
10284  |f11|d329|)|4α=↓|4|"C2,|41,|41,|41,|43|"C|'|f18|d329|)|4α=↓|
10285  4|"C1,|41,|41,|41,|41,|43|"C|'|f25|d329|)|4α=↓|4|"C1,|46,|44
10286  |"C|'>|>|f5|d329|)|4α=↓|4|"C5,|41,|44|"C|'|f12|d329|)|4α=↓|4
10290  |"c2,|42,|42,|42|"C|'|f19|d329|)|4α=↓|4|"C1,|41,|41,|49|"C|'
10292  |f26|d329|)|4α=↓|4|"C1,|48,|41,|42|"C|'>|>|f6|d329|)|4α=↓|4|
10295  "C4,|41,|45|"C|'|f13|d329|)|4α=↓|4|"C2,|44,|43|"C|'
10297  |f20|d329|)|4α=↓|4|"C1,|42,|44,|42|"C|'|f27|d329|)|4α=↓|4|"C
10298  1,|413,|42|"C|'>|>|f7|d329|)|4α=↓|4|"C4,|47|"C|'
10302  |f14|d329|)|4α=↓|4|"C2,|414|"C|'|f21|d329|)|4α=↓|4|"C1,|42,|
10303  41,|41,|41,|42|"C|'|f28|d329|)|4α=↓|4|"C1,|428|"C|'
10305  >{A6}Several things can be observed in this table.|'
10314  !|4|4|4|4a)|9As mentioned earlier, the last quotient 
10320  is always 2 or more. Furthermore, we have the 
10329  obvious identity|'{A6}|ε|"Cx|β1,|4.|4.|4.|4,|4x|βn|βα_↓|β1,|
10331  4x|βn|4α+↓|41|"C|4α=↓|4|"Cx|β1,|4.|4.|4.|4,|4x|βn|βα_↓|β1,|4
10331  x|βn,|41|"C|J!(43)|;{A6}|πwhich shows how partial 
10336  fractions whose last quotient is unity are related 
10344  to regular continued fractions.|'!|4|4|4|4b)|9The 
10349  values in the right-hand columns have a simple 
10357  relationship to the values in the left-hand columns; 
10365  can the reader see the correspondence before 
10372  reading any further? We have the identity|'{A6}|ε1|4α_↓|4|"C
10379  x|β1,|4x|β2,|4.|4.|4.|4,|4x|βn|"C|4α=↓|4|"C1,|4x|β1|4α_↓|41,
10379  |4x|β2,|4.|4.|4.|4,|4x|βn|"C;|J!(44)|;{A6}|πsee 
10381  exercise 9.|'!!c)|9There is symmetry between 
10387  left and right in the _rst two columns: If |ε|"CA|β1,|4A|β2,
10396  |4.|4.|4.|4,|4A|βt|"C |πoccurs, so does |ε|"CA|βt,|4.|4.|4.|
10400  4,|4A|β2,|4A|β1|"C. |πThis will always be the 
10406  case (see exercise 26).|'!|4|4|4|4d)|9If we examine 
10413  all of the quotients in the table, we _nd that 
10423  there are 96 in all, of which |f39|d396|)|4α=↓|440.6 
10431  percent are equal to 1, |f21|d396|)|4α=↓|421.9 
10437  percent are equal to 2, |f8|d396|)|4α=↓|48.3 
10443  percent are equal to 3; this agrees reasonably 
10451  well with the probabilities listed above.|'{A12}|≡T|≡h|≡e 
10458  |≡n|≡u|≡m|≡b|≡e|≡r |≡o|≡f |≡d|≡i|≡v|≡i|≡s|≡i|≡o|≡n 
10461  |≡s|≡t|≡e|≡p|≡s|≡.|9|4Let us now return to our 
10467  original problem and investigate |εT|βn, |πthe 
10473  average number of division steps when |εv|4α=↓|4n. 
10480  |π{H12}({H10}See Eq. (19).{H12}) {H10}Here are 
10485  some sample values of |εT|βn:|'{A6}|h|εT|βn|4|∂α=↓|4999|9999
10490  |9999|9999|99999|99999|4.|4.|4.|40000|900000|900000|E|n|;
10491  | n|4|Lα=↓|4|4|195|9|4|196|9|4|197|9|4|198|9|4|199|9100|9101
10491  |9102|9103|9104|9105>{A3}| T|βn|4|Lα=↓|45.0|94.4|95.3|94.8|9
10492  4.7|9|44.6|9|45.3|9|44.6|9|45.3|9|44.7|9|44.6>
10493  | n|4|Lα=↓|4996|9997|9998|9999|91000|91001|4|¬O|4|¬O|4|¬O|49
10493  999|910000|910001>{A3}| T|βn|4|Lα=↓|4|46.5|9|47.3|9|47.0|9|4
10494  6.8|9|9|4|16.4|9|9|4|16.7|4|¬O|4|¬O|4|¬O|4|9|48.6|9!|4|1|18.
10494  3|9!|4|1|19.1>{A3}| n|4|Lα=↓|449999|950000|950001|4|¬O|4|¬O|
10495  4|¬O|499999|9100000|9100001>{A3}| T|βn|4|Lα=↓|4|9|410.6|9!|4
10496  |1|19.7|9|9|4|110.0|4|¬O|4|¬O|4|¬O|4|9|4|110.7|9!|4|1|1|110.
10496  3|9!|4|1|1|111.0>{A6}|πNote the somewhat erratic 
10501  behavior; |εT|βn |πtends to be higher than its 
10509  neighbors when |εn |πis prime, and it is correspondingly 
10518  lower when |εn |πhas many divisors. (In this 
10526  list, 97, 101, 103, 997, and 49999 are primes; 
10535  10001|4α=↓|473|4|¬O|4137, 50001|4α=↓|43|4|¬O|47|4|¬O|42381, 
10537  99999|4α=↓|43|4|¬O|43|4|¬O|441|4|¬O|4271, 100001|4α=↓|411|4|
10538  ¬O|49091.) It is not di∃cult to understand why 
10546  this happens; if gcd(|εu,|4v)|4α=↓|4d, |πEuclid's 
10551  algorithm applied to |εu |πand |εv |πbehaves 
10558  essentially the same as if it were applied to 
10567  |εu/d |πand |εv/d. |πTherefore, when |εv|4α=↓|4n 
10573  |πhas several divisors, there are many choices 
10580  of |εu |πfor which |εn |πbehaves as if it were 
10590  smaller.|'!|4|4|4|4Accordingly let us consider 
10595  |εanother |πquantity, |ε|≤t|βn, |πwhich is the 
10601  average number of division steps when |εv|4α=↓|4n 
10608  |πand when |εu |πis |εrelatively prime |πto |εn. 
10616  |πThus|'{A6}|ε|≤t|βn|4α=↓|4|(1|d2|≤'(n)|)|4|↔k|uc|)0|¬Em|¬Wn
10617  ,|d|πgcd|ε(m,n)α=↓1|)|4T(m,|4n).|J!(45)|;{A6}|πIt 
10619  follows that|'{A6}|εT|βn|4α=↓|4|(1|d2n|)|4|↔k|uc|)d|¬Dn|)|4|
10621  ≤'(d)|≤t|βd.|J!(46)|;{A6}|πHere is a table of 
10627  |ε|≤t|βn |πfor the same values of |εn |πconsidered 
10635  above:|'{A6}|h|ε|≤t|βn|4|∂α=↓|4996|9997|9998|9999|91000|9000
10636  0|4.|4.|4.|40000|900000|900000|E|n|;| n|4|Lα=↓|4|4|195|9|4|1
10637  96|9|4|197|9|4|198|9|4|199|9100|9101|9102|9103|9104|9105>
10638  {A3}| |≤t|βn|4|Lα=↓|45.4|95.3|95.3|95.6|95.2|9|4|15.2|9|4|15
10638  .4|9|4|15.3|9|4|15.4|9|4|15.3|9|4|15.|4|¬O|4|¬O|49999|910000
10638  |910001>{A3}| |≤t|βn|4|Lα=↓|4|4|17.2|9|4|17.3|9|4|17.3|9|4|1
10639  7.3|9|9|1|4|17.3|9|9|1|4|17.4|4|¬O|4|¬O|4|¬O|4|4|19.21|9|9|1
10639  |4|19.21|9|9|1|4|19.22>{A3}| n|4|Lα=↓|449999|950000|950001|4
10640  |¬O|4|¬O|4|¬O|4|4|199999|9100000|9100001>| |≤t|βn|4|Lα=↓|4|4
10641  |110.58|9|4|110.57|9|4|110.59|4|¬O|4|¬O|4|¬O|411.170|9|4|111
10641  .172|9|4|111.172>{A6}|πClearly |ε|≤t|βn |πis 
10645  much more well-behaved than |εT|βn, |πand it 
10652  should be more susceptible to analysis. Inspection 
10659  of a table of |ε|≤t|βn |πfor small |εn |πreveals 
10668  some curious anomalies; for example, |ε|≤t|β5|β0|4α=↓|4|≤t|β
10673  1|β0|β0 |πand |ε|≤t|β6|β0|4α=↓|4|≤t|β1|β2|β0. 
10676  |πBut as |εn |πgrows, the values of |ε|≤t|βn 
10684  |πbehave quite regularly indeed, as the above 
10691  table indicates, and they show no signi_cant 
10698  relation to the factorization properties of |εn. 
10705  |πIf the reader will plot the values of |ε|≤t|βn 
10714  |πversus ln |εn |πon graph paper, for the values 
10723  of |ε|≤t|βn |πgiven above, he will see that the 
10732  values lie very nearly on a straight line, and 
10741  that the equation|'{A6}|ε|≤t|βn|4|¬V|40.843|4|πln|4|εn|4α+↓|
10744  41.47|J!(47)|;{A6}|πgives a very good approximation 
10750  to |ε|≤t|βn.|'|π!|4|4|4|4We can account for this 
10757  behavior if we return to the considerations of 
10765  Theorem W. Note that in Euclid's algorithm as 
10773  expressed in (15) we have|'{A6}|ε|(V|β0|d2U|β0|)|4|(V|β1|d2U
10778  |β1|)|4|¬O|4|¬O|4|¬O|4|(V|βt|βα_↓|β1|d2U|βt|βα_↓|β1|)|4α=↓|4
10778  |(V|βt|βα_↓|β1|d2U|β0|)|4,|;{A6}|πsince |εU|βk|βα+↓|β1|4α=↓|
10780  4V|βk; |πtherefore, if |εU|4α=↓|4U|β0 |πand |εV|4α=↓|4V|β0 
10786  |πare relatively prime, and if there are |εt 
10794  |πdivision steps, we have|'{A6}|εX|β0X|β1|4|¬O|4|¬O|4|¬O|4X|
10798  βt|βα_↓|β1|4α=↓|41/U.|;{A6}|πSetting |εU|4α=↓|4N 
10801  |πand |εV|4α=↓|4m|4|¬W|4N, |πwe _nd that|'{A6}|πln|4|εX|β0|4
10806  α+↓|4|πln|4|εX|β1|4α+↓|1|1|¬O|4|¬O|4|¬O|1|1α+↓|4|πln|4|εX|βt
10806  |βα_↓|β1|4α=↓|4α_↓|πln|4|εN.|J!(48)|;{A6}|πWe 
10808  know the approximate distribution of |εX|β0,|4X|β1,|4X|β2,|4
10813  .|4.|4.|4, |πso we can use this equation to estimate|'
10822  {A6}|εt|4α=↓|4T(N,|4m)|4α=↓|4T(m,|4N)|4α_↓|41.|;
10823  |π!|4|4|4|4Returning to the formulas preceding 
10828  Theorem W, we _nd that the average value of ln 
10838  |εX|βn, |πwhen |εX|β0 |πis a real number uniformly 
10846  distributed in [0,|41), is|'{A6}{A6}|ε{A6}|↔j|ur1|)0|)|4|πln
10850  |4|εxF|ur|↔0|)n|)(x)dx|4α=↓|4|↔j|ur1|)0|)|4|πln|4|εxf|βn(x)d
10850  x/(1|4α+↓|4x),|J!(49)|;{A12}|πwhere |εf|βn(x) 
10853  |πis de_ned in (31). Now|'{A6}|εf|βn(x)|4α=↓|4|(1|d2|πln|42|
10858  )|4α+↓|4|εO(2|gα_↓|gn),|J!(50)|;{A6}|πusing the 
10861  facts derived earlier (see exercise 23); hence 
10868  the average value of |πln|4|εX|βn |πis, to a 
10876  very good approximation,|'|Hβ*?*?*?{U0}{H10L12M29}58320#Compute